Computing the Mandelbrot Set

/**
* @file mandelbrot.cpp
* @brief Program to compute the Mandelbrot Set.
* 
* This program can compute images of the Mandelbrot Set in its
* entirety or close ups of boundary areas.
*
* get cli11.hpp at https://github.com/CLIUtils/CLI11
*
* MacOS
* brew install cli11<br>
*
* Add to .zshrc
* export CPPFLAGS="-I/opt/homebrew/include"
* export LDFLAGS="-L/opt/homebrew/lib"
*
* Compile with
* g++ $CPPFLAGS $LDFLAGS --std=gnu++20 -o mandelbrot mandelbrot.cpp
*
* Linux put CLI folder in /urs/local/include <br>
* g++ --std=gnu++20 -o mandelbrot mandelbrot.cpp
*
* @author Richard Lesh
* @date 2025-03-12
* @version 1.0
*/

#include <algorithm>
#include <CLI/CLI.hpp>
#include <cmath>
#include <cstdio>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>

const int MAX_PIXEL_LEVELS = 256;           // 256 for 8-bit

using namespace std;

// Array used to store the image as an array of arrays
// Each elemnt of the array will be a pointer to an array containing the
// pixels for a row.
unique_ptr<unique_ptr<int32_t[]>[]> image = nullptr;

/**
* @brief Compute the number of iterations to escape.
* 
* Function computes the number of iterations of the formula 
* z(n+1) = z(n)^2 + c can be executed before the value of z 
* exceeds the escape limit radius of 2.0.
* 
* @param x0 Real part of c.
* @param y0 Imaginary part of c.
* @param max_iteration Number of iterations, if exceeded, when we will declare as
*                      not escaping and therefore in the set.
* @return Value [0, max_iteration) representing the escape iterations.
*                      0 is considered as not escaping.
*/
int compute_escape(double x0, double y0, int max_iteration) {
    double x = 0;                       // We start with z(0) = 0 + 0i
    double y = 0;
    double x2 = 0;                      // Used to compute squares of the
    double y2 = 0;                      // individual components of z
    int iteration = 0;                  // Number of iterations before escaping
    
    // Loop until the magnitude |z| < 2 or we reach max_iteration
    while (x2 + y2 <= 4 && iteration < max_iteration) {
        y = 2 * x * y + y0;
        x = x2 - y2 + x0;
        x2 = x * x;
        y2 = y * y;
        ++iteration;
    }
    if (iteration == max_iteration) {   // Return 0 (black) if we are in the set
        return 0L;
    } else {
        return iteration;
    }
}

/**
* @brief Compute a portion of the Mandelbrot Set.
* 
* Function computes a portion of the Mandelbrot Set row by row.
* The row is stored in the image array.
* 
* @param row Row to compute [0, image_size)
* @param x_min Minimum real axis bound to plot.
* @param x_max Maximum real axis bound to plot.
* @param y_min Minimum imaginary axis bound to plot.
* @param y_max Maximum imaginary axis bound to plot.
* @param max_iteration Limit to the number of iterations we will compute per point.
* @param image_size Size, in pixels, of the image to generate.
*/
void compute_mandelbrot(int row, double x_min, double x_max, double y_min, double y_max, int max_iteration, int image_size) {
    double width = x_max - x_min;
    double height = y_max - y_min;
    int modulus = ceil(image_size / 40.);
    if (row < image_size) {
        image[row] = make_unique<int32_t[]>(image_size);
        double y0 = y_max - (height / image_size) * (row + double(0.5));
        for (int i = 0; i < image_size; ++i) {
            double x0 = (width / image_size) * (i + double(0.5)) + x_min;
            int level = compute_escape(x0, y0, max_iteration);
            image[row][i] = level;
        }
        if (row % modulus == 0) cout << '.' << flush;
    }
}

/**
* @brief Write out a PGM file.
* 
* Function writes out image data to a PGM 8-bit format file. This function 
* normalizes the level values using a histogram.
* 
* @param filename Pathspec to the file that will be created.
* @param max_iteration Limit to the number of iterations we will compute per point.
* @param image_size Size, in pixels, of the image to save.
*/
void output_file(string filename, int max_iteration, int image_size) {    
    vector<int> histogram(max_iteration, 0);
    string ext = ".pgm";

    cout << endl << "Creating " << (filename + ext) << "..." << endl;
    ofstream out_file((filename + ext).c_str(), ios::binary);
    if (!out_file) {                                    // Check if file opened successfully
        cerr << "Error opening file " << (filename + ext) << " for writing!" << endl;
        exit(1);
    }
    // File header is P5 for binary gray map, x_size y_size, max pixel value per RGB component
    out_file << "P5\n" << image_size << " " << image_size << "\n" << (MAX_PIXEL_LEVELS - 1) << "\n";

    // Compute histogram
    int totalCount = 0;
    for (int j = 0; j < image_size; ++j) {
        for (int i = 0; i < image_size; ++i) {
            if (image[j][i] > 0) {
                ++histogram[image[j][i]];
                ++totalCount;
            }
        }
    }
    // Compute cumulative histogram
    int count = 0;
    for (int i = 1; i < max_iteration; ++i) {
        count += histogram[i];
        histogram[i] = count;
    }
    // Normalize using histogram
    for (int j = 0; j < image_size; ++j) {
        for (int i = 0; i < image_size; ++i) {
            uint8_t level = histogram[image[j][i]] / (double) (totalCount + 1) * MAX_PIXEL_LEVELS;
            out_file.write(reinterpret_cast<const char*>(&level), 1);
        }
    }
    out_file.close();
}

/**
* @brief Main function.
* 
* The entry point of the program.
*
* Syntax: mandelbrot <x> <y> <radius> [--max_iterations=n] [--image_size=n] [--filename=name]
* 
* x & y specify the coordinates of the Mandelbrot Set to display at the center of the image.
*
* radius is the radius in Mandelbrot Set coordinates out to the extreme edges of the image.
*
* max_iterations is the number of iterations maximum per pixel to determine if it is in the 
* set.  The max_iterations will need to get larger as radius gets smaller.
*
* image_size is the size in pixels for the width and height of the image.
*
* @return int Returns 0 on successful execution.
*/
int main(int argc, char** argv) {
    CLI::App app{"Program to generate Mandelbrot images."};

    // Positional arguments: x, y, radius (as strings for now)
    vector<string> coords;
    app.add_option("coords", coords, "Location to plot (x y radius)")
        ->expected(3)  // require exactly 3 args
        ->required();

    // Optional arguments
    string filename = "image";
    int max_iteration = 10000;
    int image_size = 500;

    app.add_option("-f,--filename", filename, "Base filename (no extension)");
    app.add_option("-m,--max_iterations", max_iteration, "Max. iterations");
    app.add_option("-s,--image_size", image_size, "Image size in pixels");

    CLI11_PARSE(app, argc, argv);

    // Get the three required coordinates to plot
    double x, y, radius;
    x = stod(coords[0]);
    y = stod(coords[1]);
    radius = stod(coords[2]);

    // Check the image size
    if (image_size > sqrt(numeric_limits<int>::max())) {
        cout << "Image size can not exceed " << sqrt(numeric_limits<int>::max()) << endl;
        return 1;
    }
    
    // Compute limits of image
    double x_min = x - radius;
    double x_max = x + radius;
    double y_min = y - radius;
    double y_max = y + radius;

    // Output settings
    cout << "Maximum Iterations: " << max_iteration << endl;
    cout << "Image Bounds: " << setprecision(numeric_limits<double>::digits10)
        << x_min << ", " << x_max << ", " << y_min << ", " << y_max << endl;
    cout << "Image Size: " << image_size << " x " << image_size << endl;
    
    // Allocate an array for the row pointers.
    image = make_unique<unique_ptr<int32_t[]>[]>(image_size);
    if (image == nullptr) {                             // Check for failure
        cerr << "Memory allocation failed!" << endl;
        exit(1);
    }

    // Compute the image
    for (int i = 0; i < image_size; ++i) {
        compute_mandelbrot(i, x_min, x_max, y_min, y_max, max_iteration, image_size);
    }

    // Generate the output file
    output_file(filename, max_iteration, image_size);
    return 0;
}

/**
* @file mandelbrot_threaded.cpp
* @brief Program to compute the Mandelbrot Set.
* 
* This program can compute images of the Mandelbrot Set in its
* entirety or close ups of boundary areas.  It will use as many
* CPUs as are available to speed up the rendering.
*
* get cli11.hpp at https://github.com/CLIUtils/CLI11
*
* MacOS
* brew install cli11<br>
*
* Add to .zshrc
* export CPPFLAGS="-I/opt/homebrew/include"
* export LDFLAGS="-L/opt/homebrew/lib"
*
* Compile with
* g++ $CPPFLAGS $LDFLAGS --std=gnu++20 -o mandelbrot_threaded mandelbrot_threaded.cpp
*
* Linux put CLI folder in /urs/local/include <br>
* g++ --std=gnu++20 -o mandelbrot_threaded mandelbrot_threaded.cpp
*
* @author Richard Lesh
* @date 2025-03-12
* @version 1.0
*/

#include <algorithm>
#include <CLI/CLI.hpp>
#include <cmath>
#include <cstdio>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <string>
#include <thread>
#include <vector>

const int MAX_PIXEL_LEVELS = 256;           // 256 for 8-bit

using namespace std;

// Array used to store the image in a thread safe way as an array of arrays
// Each elemnt of the array will be a pointer to an array containing the
// pixels for a row.  Rows will be computed by threads.
unique_ptr<unique_ptr<int32_t[]>[]> image = nullptr;
// Atomic int used to assign tasks (row to compute) to the threads
atomic_int counter(0);

/**
* @brief Compute the number of iterations to escape.
* 
* Function computes the number of iterations of the formula 
* z(n+1) = z(n)^2 + c can be executed before the value of z 
* exceeds the escape limit radius of 2.0.
* 
* @param x0 Real part of c.
* @param y0 Imaginary part of c.
* @param max_iteration Number of iterations, if exceeded, when we will declare as
*                      not escaping and therefore in the set.
* @return Value [0, max_iteration) representing the escape iterations.
*                      0 is considered as not escaping.
*/
int compute_escape(double x0, double y0, int max_iteration) {
    double x = 0;                       // We start with z(0) = 0 + 0i
    double y = 0;
    double x2 = 0;                      // Used to compute squares of the
    double y2 = 0;                      // individual components of z
    int iteration = 0;                  // Number of iterations before escaping
    
    // Loop until the magnitude |z| < 2 or we reach max_iteration
    while (x2 + y2 <= 4 && iteration < max_iteration) {
        y = 2 * x * y + y0;
        x = x2 - y2 + x0;
        x2 = x * x;
        y2 = y * y;
        ++iteration;
    }
    if (iteration == max_iteration) {   // Return 0 (black) if we are in the set
        return 0L;
    } else {
        return iteration;
    }
}

/**
* @brief Compute a portion of the Mandelbrot Set.
* 
* Function computes a portion of the Mandelbrot Set row by row.
* The atomic integer counter is used to get the row to compute.
* The row is stored in the image array.
* 
* @param x_min Minimum real axis bound to plot.
* @param x_max Maximum real axis bound to plot.
* @param y_min Minimum imaginary axis bound to plot.
* @param y_max Maximum imaginary axis bound to plot.
* @param max_iteration Limit to the number of iterations we will compute per point.
* @param image_size Size, in pixels, of the image to generate.
*/
void compute_mandelbrot(double x_min, double x_max, double y_min, double y_max, int max_iteration, int image_size) {
    double width = x_max - x_min;
    double height = y_max - y_min;
    int modulus = ceil(image_size / 40.);
    while (true) {
        int j = counter.fetch_add(1);
        if (j < image_size) {
            image[j] = make_unique<int32_t[]>(image_size);
            double y0 = y_max - (height / image_size) * (j + double(0.5));
            for (int i = 0; i < image_size; ++i) {
                double x0 = (width / image_size) * (i + double(0.5)) + x_min;
                int level = compute_escape(x0, y0, max_iteration);
                image[j][i] = level;
            }
            if (j % modulus == 0) cout << '.' << flush;
        } else {
            break;
        }
    }
}

/**
* @brief Write out a PGM file.
* 
* Function writes out image data to a PGM 8-bit format file. This function 
* normalizes the level values using a histogram.
* 
* @param filename Pathspec to the file that will be created.
* @param max_iteration Limit to the number of iterations we will compute per point.
* @param image_size Size, in pixels, of the image to save.
*/
void output_file(string filename, int max_iteration, int image_size) {    
    vector<int> histogram(max_iteration, 0);
    string ext = ".pgm";

    cout << endl << "Creating " << (filename + ext) << "..." << endl;
    ofstream out_file((filename + ext).c_str(), ios::binary);
    if (!out_file) {                                    // Check if file opened successfully
        cerr << "Error opening file " << (filename + ext) << " for writing!" << endl;
        exit(1);
    }
    // File header is P5 for binary gray map, x_size y_size, max pixel value per RGB component
    out_file << "P5\n" << image_size << " " << image_size << "\n" << (MAX_PIXEL_LEVELS - 1) << "\n";

    // Compute histogram
    int totalCount = 0;
    for (int j = 0; j < image_size; ++j) {
        for (int i = 0; i < image_size; ++i) {
            if (image[j][i] > 0) {
                ++histogram[image[j][i]];
                ++totalCount;
            }
        }
    }
    // Compute cumulative histogram
    int count = 0;
    for (int i = 1; i < max_iteration; ++i) {
        count += histogram[i];
        histogram[i] = count;
    }
    // Normalize using histogram
    for (int j = 0; j < image_size; ++j) {
        for (int i = 0; i < image_size; ++i) {
            uint8_t level = histogram[image[j][i]] / (double) (totalCount + 1) * MAX_PIXEL_LEVELS;
            out_file.write(reinterpret_cast<const char*>(&level), 1);
        }
    }
    out_file.close();
}

/**
* @brief Main function.
* 
* The entry point of the program.
*
* Syntax: mandelbrot <x> <y> <radius> [--max_iterations=n] [--image_size=n] [--filename=name]
* 
* x & y specify the coordinates of the Mandelbrot Set to display at the center of the image.
*
* radius is the radius in Mandelbrot Set coordinates out to the extreme edges of the image.
*
* max_iterations is the number of iterations maximum per pixel to determine if it is in the 
* set.  The max_iterations will need to get larger as radius gets smaller.
*
* image_size is the size in pixels for the width and height of the image.
*
* @return int Returns 0 on successful execution.
*/
int main(int argc, char** argv) {
    CLI::App app{"Program to generate Mandelbrot images."};
    // use cxxpots to process command line options.
    vector<string> coords;
    app.add_option("coords", coords, "Location to plot (x y radius)")
        ->expected(3)  // require exactly 3 args
        ->required();

    // Optional arguments
    string filename = "image";
    int max_iteration = 10000;
    int image_size = 500;

    app.add_option("-f,--filename", filename, "Base filename (no extension)");
    app.add_option("-m,--max_iterations", max_iteration, "Max. iterations");
    app.add_option("-s,--image_size", image_size, "Image size in pixels");

    CLI11_PARSE(app, argc, argv);

    // Get the three required coordinates to plot
    double x, y, radius;
    x = stod(coords[0]);
    y = stod(coords[1]);
    radius = stod(coords[2]);

    // Check the image size
    if (image_size > sqrt(numeric_limits<int>::max())) {
        cout << "Image size can not exceed " << sqrt(numeric_limits<int>::max()) << endl;
        return 1;
    }
    
    // Compute limits of image
    double x_min = x - radius;
    double x_max = x + radius;
    double y_min = y - radius;
    double y_max = y + radius;

    // Output settings
    cout << "Maximum Iterations: " << max_iteration << endl;
    cout << "Image Bounds: " << setprecision(numeric_limits<double>::digits10)
        << x_min << ", " << x_max << ", " << y_min << ", " << y_max << endl;
    cout << "Image Size: " << image_size << " x " << image_size << endl;
    
    // Allocate an array for the row pointers.
    image = make_unique<unique_ptr<int32_t[]>[]>(image_size);
    if (image == nullptr) {                             // Check for failure
        cerr << "Memory allocation failed!" << endl;
        exit(1);
    }

    // Get the number of CPU cores.
    unsigned int num_cpus = thread::hardware_concurrency();
    if (num_cpus == 0) {
        cout << "Could not determine number of CPUs.  Using 4.\n";
        num_cpus = 4;
    } else {
        cout << "Number of CPU cores: " << num_cpus << endl;
    }

    // Create the threads
    vector<thread> threads;
    for (int i = 0; i < num_cpus; ++i) {
        threads.emplace_back(compute_mandelbrot, x_min, x_max, y_min, y_max, max_iteration, image_size);
    }
    // Wait for all threads to complete
    for (auto& t : threads) {
        t.join();  
    }

    // Generate the output file
    output_file(filename, max_iteration, image_size);
    return 0;
}

/**
* @file mandelbrot_threaded_color.cpp
* @brief Program to compute the Mandelbrot Set.
* 
* This program can compute images of the Mandelbrot Set in its
* entirety or close ups of boundary areas.  It will use as many
* CPUs as are available to speed up the rendering.
*
* It also allows two RGB colors to be specified to define the
* coloring of the set based on the escape value of a pixel.
* Colors can be specified using #xxx or #xxxxxx hexidecimal.
*
* get cli11.hpp at https://github.com/CLIUtils/CLI11
*
* MacOS:
* brew install cli11
*
* Add to .zshrc<br>
* export CPPFLAGS="-I/opt/homebrew/include"<br>
* export LDFLAGS="-L/opt/homebrew/lib"
*
* Compile with<br>
* g++ $CPPFLAGS $LDFLAGS --std=gnu++20 -o mandelbrot_threaded_color mandelbrot_threaded_color.cpp
*
* Linux:<br>
* put CLI folder in /urs/local/include <br>
* g++ --std=gnu++20 -o mandelbrot_threaded_color mandelbrot_threaded_color.cpp
*
* @author Richard Lesh
* @date 2025-03-12
* @version 1.0
*/

#include <algorithm>
#include <CLI/CLI.hpp>
#include <cmath>
#include <cstdio>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <string>
#include <thread>
#include <vector>

const int MAX_PIXEL_LEVELS = 256;           // 256 or (256 * 256) for 8-bit or 16-bit respectively
bool USE_COLOR = true;                      // True to allow colors to be used
#undef USE_RAINBOW                          // Define to use rainbow colors

#if MAX_PIXEL_LEVELS > 256
    #define PIXEL_TYPE uint16_t
#else
    #define PIXEL_TYPE uint8_t
#endif

using namespace std;

// Array used to store the image in a thread safe way as an array of arrays
// Each elemnt of the array will be a pointer to an array containing the
// pixels for a row.  Rows will be computed by threads.
unique_ptr<unique_ptr<int32_t[]>[]> image = nullptr;
// Atomic int used to assign tasks (row to compute) to the threads
atomic_int counter(0);

/**
* @brief Stucture for storing RGB colors.
* 
* Stores the RGB components of an RGB color.
* Can be 8-bit or 16-bit per channel.
*/
template <typename Depth>
struct Color {
    Depth R;
    Depth G;
    Depth B;

    /**
    * @brief Default Color constructor.
    * 
    * Creates a Color object set to white.
    */
    Color() {
        R = G = B = MAX_PIXEL_LEVELS - 1;
    }
    
    /**
    * @brief Constructor for RGB colors.
    * 
    * Each color compontent can be in the range [0, MAX_PIXEL_LEVELS)
    * @param r Red component
    * @param g Green component
    * @param b Blue component
    */
    Color(int r, int g, int b) {
        R = max(0, min(r, MAX_PIXEL_LEVELS - 1));
        G = max(0, min(g, MAX_PIXEL_LEVELS - 1));
        B = max(0, min(b, MAX_PIXEL_LEVELS - 1));
    }

    /**
    * @brief Constructor using hexadecimal string.
    * 
    * Each color compontent is specified in hexadecimal in string form
    * "#xxx"  or "#xxxxxx".
    *
    * @param hex Hexadecimal color string
    */
    Color(string hex) {
        // Ensure the input is valid
        string hex_color = hex;
        if (hex_color[0] == '#') hex_color = hex_color.substr(1);   // Remove `#` if present

        if (hex_color.length() == 3) {                              // Shorthand #RGB format
            hex_color = {hex_color[0], hex_color[0], hex_color[1], hex_color[1], hex_color[2], hex_color[2]};
        }

        if (hex_color.length() != 6) {
            throw invalid_argument("Invalid hex color format");
        }

        // Convert hex to RGB values
        R = static_cast<Depth>(stoi(hex_color.substr(0, 2), nullptr, 16));
        G = static_cast<Depth>(stoi(hex_color.substr(2, 2), nullptr, 16));
        B = static_cast<Depth>(stoi(hex_color.substr(4, 2), nullptr, 16));
    }
    bool operator==(const Color<Depth> &arg) const {
        return R == arg.R && G == arg.G && B == arg.B;
    }
    bool isGray() {
        return R == G && G == B;
    }
    friend ostream;
};

template <typename Depth>
ostream &operator<<(ostream &os, const Color<Depth> &rgb) {
    os << "RGB(" << (int)rgb.R << "," << (int)rgb.G << "," << (int)rgb.B << ")";
    return os;
}

typedef Color<PIXEL_TYPE> RGB;

const RGB BLACK{0,0,0};
const RGB WHITE{};

/**
* @brief Function to blend two RGB colors.
* 
* Function to blend two RGB colors based on a level in [0.0, 1.0)
* Exactly 0.0 will return black.
*
* @param level Level value to control blending [0.0, 1.0)
* @param color1 Color for near 0.0 level
* @param color2 Color for near 1.0 level
* return Blended color
*/
RGB blend_colors(double level, const RGB& color1, const RGB& color2) {
    if (level <= 0.0) return BLACK;
    if (level >= 1.0) return color2;

    return RGB(
        static_cast<unsigned int>(color1.R + level * (color2.R - color1.R + 1)),
        static_cast<unsigned int>(color1.G + level * (color2.G - color1.G + 1)),
        static_cast<unsigned int>(color1.B + level * (color2.B - color1.B + 1))
    );
}

/**
* @brief Function to generate rainbow colors.
* 
* Function to a rainbow of colors based on a level in [0.0, 1.0)
* Exactly 0.0 will return black.
*
* @param level Level value to control rainbow [0.0, 1.0)
* return Rainbow color
*/
RGB level_2_RGB(double level) {
    if (level <= 0.0) return BLACK;
    double H = (log(level + 1)/log(log(2.0)) - 1);
    double R = (sin(2.0 * M_PI * H) + 1.0) / 2.0;
    double G = (sin(2.0 * M_PI * H + 2.0 * M_PI / 3.0) + 1.0) / 2.0;
    double B = (sin(2.0 * M_PI * H + 4.0 * M_PI / 3.0) + 1.0) / 2.0;
    RGB result((int)floor(R * MAX_PIXEL_LEVELS), (int)floor(G * MAX_PIXEL_LEVELS), (int)floor(B * MAX_PIXEL_LEVELS));
    return result;
}

/**
* @brief Compute the number of iterations to escape.
* 
* Function computes the number of iterations of the formula 
* z(n+1) = z(n)^2 + c can be executed before the value of z 
* exceeds the escape limit radius of 2.0.
* 
* @param x0 Real part of c.
* @param y0 Imaginary part of c.
* @param max_iteration Number of iterations, if exceeded, when we will declare as
*                      not escaping and therefore in the set.
* @return Value [0, max_iteration) representing the escape iterations.
*                      0 is considered as not escaping.
*/
int compute_escape(double x0, double y0, int max_iteration) {
    double x = 0;                       // We start with z(0) = 0 + 0i
    double y = 0;
    double x2 = 0;                      // Used to compute squares of the
    double y2 = 0;                      // individual components of z
    int iteration = 0;                  // Number of iterations before escaping
    
    // Loop until the magnitude |z| < 2 or we reach max_iteration
    while (x2 + y2 <= 4 && iteration < max_iteration) {
        y = 2 * x * y + y0;
        x = x2 - y2 + x0;
        x2 = x * x;
        y2 = y * y;
        ++iteration;
    }
    if (iteration == max_iteration) {   // Return 0 (black) if we are in the set
        return 0L;
    } else {
        return iteration;
    }
}

/**
* @brief Compute a portion of the Mandelbrot Set.
* 
* Function computes a portion of the Mandelbrot Set row by row.
* The atomic integer counter is used to get the row to compute.
* The row is stored in the image array.
* 
* @param x_min Minimum real axis bound to plot.
* @param x_max Maximum real axis bound to plot.
* @param y_min Minimum imaginary axis bound to plot.
* @param y_max Maximum imaginary axis bound to plot.
* @param max_iteration Limit to the number of iterations we will compute per point.
* @param image_size Size, in pixels, of the image to generate.
*/
void compute_mandelbrot(double x_min, double x_max, double y_min, double y_max, int max_iteration, int image_size) {
    double width = x_max - x_min;
    double height = y_max - y_min;
    int modulus = ceil(image_size / 40.);
    while (true) {
        int j = counter.fetch_add(1);
        if (j < image_size) {
            image[j] = make_unique<int32_t[]>(image_size);
            double y0 = y_max - (height / image_size) * (j + double(0.5));
            for (int i = 0; i < image_size; ++i) {
                double x0 = (width / image_size) * (i + double(0.5)) + x_min;
                int level = compute_escape(x0, y0, max_iteration);
                image[j][i] = level;
            }
            if (j % modulus == 0) cout << '.' << flush;
        } else {
            break;
        }
    }
}

/**
* @brief Write out a PGM/PPM file.
* 
* Function writes out image data to a PGM or PPM format file.  Can be an 
* 8-bit/16-bit PGM or an 8-bit/16-bit PPM.  This function normalizes the
* level values using a histogram.
* 
* @param filename Pathspec to the file that will be created.
* @param max_iteration Limit to the number of iterations we will compute per point.
* @param image_size Size, in pixels, of the image to save.
* @param color1 First color to use for near 1 iteration before escaping.
* @param color2 Second color to use for near max_iterations iterations before escaping.
*/
void output_file(string filename, int max_iteration, int image_size, RGB color1, RGB color2) {    
    vector<int> histogram(max_iteration, 0);
    string ext = ".pgm";
    if (USE_COLOR) ext = ".ppm";

    cout << endl << "Creating " << (filename + ext) << "..." << endl;
    ofstream out_file((filename + ext).c_str(), ios::binary);
    if (!out_file) {                                    // Check if file opened successfully
        cerr << "Error opening file " << (filename + ext) << " for writing!" << endl;
        exit(1);
    }
    // File header is P6 for binary pix map, x_size y_size, max pixel value per RGB component
    if (USE_COLOR) {
        out_file << "P6\n" << image_size << " " << image_size << "\n" << (MAX_PIXEL_LEVELS - 1) << "\n";
    } else {
        out_file << "P5\n" << image_size << " " << image_size << "\n" << (MAX_PIXEL_LEVELS - 1) << "\n";
    }

    // Compute histogram
    int totalCount = 0;
    for (int j = 0; j < image_size; ++j) {
        for (int i = 0; i < image_size; ++i) {
            if (image[j][i] > 0) {
                ++histogram[image[j][i]];
                ++totalCount;
            }
        }
    }
    // Compute cumulative histogram
    int count = 0;
    for (int i = 1; i < max_iteration; ++i) {
        count += histogram[i];
        histogram[i] = count;
    }
    // Normalize using histogram
    for (int j = 0; j < image_size; ++j) {
        for (int i = 0; i < image_size; ++i) {
            if (USE_COLOR) {
                RGB rgb = BLACK;
                if (image[j][i] != 0) {
#ifdef USE_RAINBOW
                    rgb = level_2_RGB(histogram[image[j][i]] / (double) (totalCount + 1));
#else
                    rgb = blend_colors(histogram[image[j][i]] / (double) (totalCount + 1), color1, color2);
#endif
                }
#if MAX_PIXEL_LEVELS > 256
                out_file.write(reinterpret_cast<const char*>(&rgb.R) + 1, 1);
                out_file.write(reinterpret_cast<const char*>(&rgb.R), 1);
                out_file.write(reinterpret_cast<const char*>(&rgb.G) + 1, 1);
                out_file.write(reinterpret_cast<const char*>(&rgb.G), 1);
                out_file.write(reinterpret_cast<const char*>(&rgb.B) + 1, 1);
                out_file.write(reinterpret_cast<const char*>(&rgb.B), 1);
#else
                out_file.write(reinterpret_cast<const char*>(&rgb.R), 1);
                out_file.write(reinterpret_cast<const char*>(&rgb.G), 1);
                out_file.write(reinterpret_cast<const char*>(&rgb.B), 1);
#endif
            } else {
#if MAX_PIXEL_LEVELS > 256
                uint16_t level = 0;
                if (image[j][i] != 0)
                    level = color1.R + histogram[image[j][i]] / (double) (totalCount + 1) * (color2.R - color1.R + 1);

                out_file.write(reinterpret_cast<const char*>(&level) + 1, 1);
                out_file.write(reinterpret_cast<const char*>(&level), 1);
#else
                uint8_t level = 0;
                if (image[j][i] != 0)
                    level = color1.R + histogram[image[j][i]] / (double) (totalCount + 1) * (color2.R - color1.R + 1);
                out_file.write(reinterpret_cast<const char*>(&level), 1);
#endif
            }
        }
    }
    out_file.close();
}

/**
* @brief Main function.
* 
* The entry point of the program.
*
* Syntax: mandelbrot <x> <y> <radius> [--max_iterations=n] [--image_size=n] 
*   [--filename=name] [--high_color=color] [--low_color=color]
*
* x & y specify the coordinates of the Mandelbrot Set to display at the center of the image.
*
* radius is the radius in Mandelbrot Set coordinates out to the extreme edges of the image.
*
* max_iterations is the number of iterations maximum per pixel to determine if it is in the 
* set.  The max_iterations will need to get larger as radius gets smaller.
*
* image_size is the size in pixels for the width and height of the image.
*
* high_color and low_color are hexadecimal color values.
*
* filename is the PGM or PPM file that will be created with the image.
*
* @return int Returns 0 on successful execution.
*/
int main(int argc, char** argv) {
    CLI::App app{"Program to generate Mandelbrot images."};
    app.set_help_flag("", "");
    // use CLI11 to process command line options.
    vector<string> coords;
    app.add_option("coords", coords, "Location to plot (x y radius)")
        ->expected(3);

    // Optional arguments
    string filename = "image";
    int max_iteration = 10000;
    int image_size = 500;
    string high_color = "#FFF";
    string low_color = "#000";
    bool show_help = false;

    app.add_flag("--help", show_help, "Show help");
    app.add_option("-f,--filename", filename, "Base filename (no extension)");
    app.add_option("-m,--max_iterations", max_iteration, "Max. iterations");
    app.add_option("-s,--image_size", image_size, "Image size in pixels");
    app.add_option("-h,--high_color", high_color, "High color hex code (#xxx or #xxxxxx)");
    app.add_option("-l,--low_color", low_color, "Low color hex code (#xxx or #xxxxxx)");

    CLI11_PARSE(app, argc, argv);

    if (show_help) {
        cout << app.help() << endl;
        return 0;
    }

    // Get the three required coordinates to plot
    double x, y, radius;
    if (coords.size() == 3) {
        x = stod(coords[0]);
        y = stod(coords[1]);
        radius = stod(coords[2]);
    } else {
        cout << "Arguments x, y and radius required but not found!" << endl;
        return 1;
    }

    // Check the image size
    if (image_size > sqrt(numeric_limits<int>::max())) {
        cout << "Image size can not exceed " << sqrt(numeric_limits<int>::max()) << endl;
        return 1;
    }
    
    // Get the low and high colors
    RGB high_RGB = RGB(high_color);
    RGB low_RGB = RGB(low_color);
    if (high_RGB.isGray() && low_RGB.isGray()) USE_COLOR = false;

    // Compute limits of image
    double x_min = x - radius;
    double x_max = x + radius;
    double y_min = y - radius;
    double y_max = y + radius;

    // Output settings
    cout << "Maximum Iterations: " << max_iteration << endl;
    cout << "High escape color: " << high_color << endl;
    cout << "Low escape color: " << low_color << endl;
    cout << (USE_COLOR ? "Using color" : "Using grayscale") << endl;
    cout << "Image Bounds: " << setprecision(numeric_limits<double>::digits10)
        << x_min << ", " << x_max << ", " << y_min << ", " << y_max << endl;
    cout << "Image Size: " << image_size << " x " << image_size << endl;
    
    // Allocate an array for the row pointers.
    image = make_unique<unique_ptr<int32_t[]>[]>(image_size);
    if (image == nullptr) {                             // Check for failure
        cerr << "Memory allocation failed!" << endl;
        exit(1);
    }

    // Get the number of CPU cores.
    unsigned int num_cpus = thread::hardware_concurrency();
    if (num_cpus == 0) {
        cout << "Could not determine number of CPUs.  Using 4.\n";
        num_cpus = 4;
    } else {
        cout << "Number of CPU cores: " << num_cpus << endl;
    }

    // Create the threads
    vector<thread> threads;
    for (int i = 0; i < num_cpus; ++i) {
        threads.emplace_back(compute_mandelbrot, x_min, x_max, y_min, y_max, max_iteration, image_size);
    }
    // Wait for all threads to complete
    for (auto& t : threads) {
        t.join();  
    }

    // Generate the output file
    output_file(filename, max_iteration, image_size, low_RGB, high_RGB);
    return 0;
}

/**
* @file mandelbrot_threaded_color_mp.cpp
* @brief Program to compute the Mandelbrot Set.
* 
* This program can compute images of the Mandelbrot Set in its
* entirety or close ups of boundary areas.  It will use as many
* CPUs as are available to speed up the rendering.  This program
* will use extended precision floating point.  The Gnu Multi-Precision
* library is required if hardware extended precision FP is not
* available.
*
* It also allows two RGB colors to be specified to define the
* coloring of the set based on the escape value of a pixel.
* Colors can be specified using #xxx or #xxxxxx hexidecimal.
*
* If 128-bit doubles are available on your platform <br>
* compile with command: g++ --std=gnu++20 -lquadmath mandelbrot_threaded_color_mp.cpp -o mandelbrot_threaded_color_mp
*
* To compile with Gnu Multi-precision floats
*
* get mpreal.h at https://github.com/advanpix/mpreal.git
*
* get cli11.hpp at https://github.com/CLIUtils/CLI11
*
* MacOS:
* put mpreal.h in /opt/homebrew/include <br>
* brew install gmp <br>
* brew install mpfr <br>
* brew install cli11
*
* Add to .zshrc<br>
* export CPPFLAGS="-I/opt/homebrew/include"<br>
* export LDFLAGS="-L/opt/homebrew/lib"
*
* Compile with<br>
* g++ $CPPFLAGS $LDFLAGS --std=gnu++20 -o mandelbrot_threaded_color_mp mandelbrot_threaded_color_mp.cpp -lmpfr -lgmp
*
* Linux:<br>
* put mpreal.h in /urs/local/include <br>
* put CLI folder in /urs/local/include <br>
* sudo apt install libgmp-dev <br>
* sudo apt install libmpfr-dev <br>
* g++ --std=gnu++20 -o mandelbrot_threaded_color_mp mandelbrot_threaded_color_mp.cpp -lmpfr -lgmp
*
* Test multi-precision with command <br>
* ./mandelbrot_threaded_color 0.3663629834227643413289327308 0.59153377326144522757749349375 1e-27 -m 100000 -d 30
*
* @author Richard Lesh
* @date 2025-03-12
* @version 1.0
*/

#include <algorithm>
#include <CLI/CLI.hpp>
#include <cmath>
#include <cstdio>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <string>
#include <thread>
#include <vector>

const int MAX_PIXEL_LEVELS = 256;           // 256 or (256 * 256) for 8-bit or 16-bit respectively
bool USE_COLOR = true;                      // True to allow colors to be used
#undef USE_RAINBOW                          // Define to use rainbow colors
#define use_extended_fp                     // Define to use extended precision floating point

#if MAX_PIXEL_LEVELS > 256
    #define PIXEL_TYPE uint16_t
#else
    #define PIXEL_TYPE uint8_t
#endif

int DIGITS = 16;                                    // Decimal digits of accuracy
#ifdef use_extended_fp
    #ifdef __FLOAT128__
        typedef __float128 double_type;
    #elif defined(__SIZEOF_LONG_DOUBLE__) && (__SIZEOF_LONG_DOUBLE__ > 8)
        #define LONG_DOUBLE
        typedef long double double_type;
    #else
        #define MPREAL
        #include <mpreal.h>                         // MPFR C++ Wrapper
        using namespace mpfr;
        #define PRECISION int(ceil(DIGITS * 3.3))   // Precision in bits (~3.3 bits per decimal digit)
        #define DOUBLE_DIGITS DIGITS
        typedef mpreal double_type;
    #endif
#else
    typedef double double_type;
#endif
#ifndef MPREAL
    #define DOUBLE_DIGITS numeric_limits<double_type>::digits10
    #define PRECISION numeric_limits<double_type>::digits
#endif

using namespace std;

// Array used to store the image in a thread safe way as an array of arrays
// Each elemnt of the array will be a pointer to an array containing the
// pixels for a row.  Rows will be computed by threads.
unique_ptr<unique_ptr<int32_t[]>[]> image = nullptr;
// Atomic int used to assign tasks (row to compute) to the threads
atomic_int counter(0);

/**
* @brief Stucture for storing RGB colors.
* 
* Stores the RGB components of an RGB color.
* Can be 8-bit or 16-bit per channel.
*/
template <typename Depth>
struct Color {
    Depth R;
    Depth G;
    Depth B;

    /**
    * @brief Default Color constructor.
    * 
    * Creates a Color object set to white.
    */
    Color() {
        R = G = B = MAX_PIXEL_LEVELS - 1;
    }
    
    /**
    * @brief Constructor for RGB colors.
    * 
    * Each color compontent can be in the range [0, MAX_PIXEL_LEVELS)
    * @param r Red component
    * @param g Green component
    * @param b Blue component
    */
    Color(int r, int g, int b) {
        R = max(0, min(r, MAX_PIXEL_LEVELS - 1));
        G = max(0, min(g, MAX_PIXEL_LEVELS - 1));
        B = max(0, min(b, MAX_PIXEL_LEVELS - 1));
    }

    /**
    * @brief Constructor using hexadecimal string.
    * 
    * Each color compontent is specified in hexadecimal in string form
    * "#xxx"  or "#xxxxxx".
    *
    * @param hex Hexadecimal color string
    */
    Color(string hex) {
        // Ensure the input is valid
        string hex_color = hex;
        if (hex_color[0] == '#') hex_color = hex_color.substr(1);   // Remove `#` if present

        if (hex_color.length() == 3) {                              // Shorthand #RGB format
            hex_color = {hex_color[0], hex_color[0], hex_color[1], hex_color[1], hex_color[2], hex_color[2]};
        }

        if (hex_color.length() != 6) {
            throw invalid_argument("Invalid hex color format");
        }

        // Convert hex to RGB values
        R = static_cast<Depth>(stoi(hex_color.substr(0, 2), nullptr, 16));
        G = static_cast<Depth>(stoi(hex_color.substr(2, 2), nullptr, 16));
        B = static_cast<Depth>(stoi(hex_color.substr(4, 2), nullptr, 16));
    }
    bool operator==(const Color<Depth> &arg) const {
        return R == arg.R && G == arg.G && B == arg.B;
    }
    bool isGray() {
        return R == G && G == B;
    }
    friend ostream;
};

template <typename Depth>
ostream &operator<<(ostream &os, const Color<Depth> &rgb) {
    os << "RGB(" << (int)rgb.R << "," << (int)rgb.G << "," << (int)rgb.B << ")";
    return os;
}

typedef Color<PIXEL_TYPE> RGB;

const RGB BLACK{0,0,0};
const RGB WHITE{};

/**
* @brief Function to blend two RGB colors.
* 
* Function to blend two RGB colors based on a level in [0.0, 1.0)
* Exactly 0.0 will return black.
*
* @param level Level value to control blending [0.0, 1.0)
* @param color1 Color for near 0.0 level
* @param color2 Color for near 1.0 level
* return Blended color
*/
RGB blend_colors(double level, const RGB& color1, const RGB& color2) {
    if (level <= 0.0) return BLACK;
    if (level >= 1.0) return color2;

    return RGB(
        static_cast<unsigned int>(color1.R + level * (color2.R - color1.R + 1)),
        static_cast<unsigned int>(color1.G + level * (color2.G - color1.G + 1)),
        static_cast<unsigned int>(color1.B + level * (color2.B - color1.B + 1))
    );
}

/**
* @brief Function to generate rainbow colors.
* 
* Function to a rainbow of colors based on a level in [0.0, 1.0)
* Exactly 0.0 will return black.
*
* @param level Level value to control rainbow [0.0, 1.0)
* return Rainbow color
*/
RGB level_2_RGB(double level) {
    if (level <= 0.0) return BLACK;
    double H = (log(level + 1)/log(log(2.0)) - 1);
    double R = (sin(2.0 * M_PI * H) + 1.0) / 2.0;
    double G = (sin(2.0 * M_PI * H + 2.0 * M_PI / 3.0) + 1.0) / 2.0;
    double B = (sin(2.0 * M_PI * H + 4.0 * M_PI / 3.0) + 1.0) / 2.0;
    RGB result((int)floor(R * MAX_PIXEL_LEVELS), (int)floor(G * MAX_PIXEL_LEVELS), (int)floor(B * MAX_PIXEL_LEVELS));
    return result;
}

/**
* @brief Compute the number of iterations to escape.
* 
* Function computes the number of iterations of the formula 
* z(n+1) = z(n)^2 + c can be executed before the value of z 
* exceeds the escape limit radius of 2.0.
* 
* @param x0 Real part of c.
* @param y0 Imaginary part of c.
* @param max_iteration Number of iterations, if exceeded, when we will declare as
*                      not escaping and therefore in the set.
* @return Value [0, max_iteration) representing the escape iterations.
*                      0 is considered as not escaping.
*/
int compute_escape(double_type x0, double_type y0, int max_iteration) {
    double_type x = 0;                      // We start with z(0) = 0 + 0i
    double_type y = 0;
    double_type x2 = 0;                     // Used to compute squares of the
    double_type y2 = 0;                     // individual components of z
    int iteration = 0;                      // Number of iterations before escaping
    
    // Loop until the magnitude |z| < 2 or we reach max_iteration
    while (x2 + y2 <= 4 && iteration < max_iteration) {
        y = 2 * x * y + y0;
        x = x2 - y2 + x0;
        x2 = x * x;
        y2 = y * y;
        ++iteration;
    }
    if (iteration == max_iteration) {       // Return 0 (black) if we are in the set
        return 0L;
    } else {
        return iteration;
    }
}

/**
* @brief Compute a portion of the Mandelbrot Set.
* 
* Function computes a portion of the Mandelbrot Set row by row.
* The atomic integer counter is used to get the row to compute.
* The row is stored in the image array.
* 
* @param x_min Minimum real axis bound to plot.
* @param x_max Maximum real axis bound to plot.
* @param y_min Minimum imaginary axis bound to plot.
* @param y_max Maximum imaginary axis bound to plot.
* @param max_iteration Limit to the number of iterations we will compute per point.
* @param image_size Size, in pixels, of the image to generate.
*/
void compute_mandelbrot(double_type x_min, double_type x_max, double_type y_min, double_type y_max, int max_iteration, int image_size) {
#ifdef MPREAL
    mpreal::set_default_prec(PRECISION);      // Set default precision in bits for each thread
#endif
    double_type width = x_max - x_min;
    double_type height = y_max - y_min;
    int modulus = ceil(image_size / 40.);
    while (true) {
        int j = counter.fetch_add(1);
        if (j < image_size) {
            image[j] = make_unique<int32_t[]>(image_size);
            double_type y0 = y_max - (height / image_size) * (j + double_type(0.5));
            for (int i = 0; i < image_size; ++i) {
                double_type x0 = (width / image_size) * (i + double_type(0.5)) + x_min;
                int level = compute_escape(x0, y0, max_iteration);
                image[j][i] = level;
            }
            if (j % modulus == 0) cout << '.' << flush;
        } else {
            break;
        }
    }
}

/**
* @brief Write out a PGM/PPM file.
* 
* Function writes out image data to a PGM or PPM format file.  Can be an 
* 8-bit/16-bit PGM or an 8-bit/16-bit PPM.  This function normalizes the
* level values using a histogram.
* 
* @param filename Pathspec to the file that will be created.
* @param max_iteration Limit to the number of iterations we will compute per point.
* @param image_size Size, in pixels, of the image to save.
* @param color1 First color to use for near 1 iteration before escaping.
* @param color2 Second color to use for near max_iterations iterations before escaping.
*/
void output_file(string filename, int max_iteration, int image_size, RGB color1, RGB color2) {    
    vector<int> histogram(max_iteration, 0);
    string ext = ".pgm";
    if (USE_COLOR) ext = ".ppm";

    cout << endl << "Creating " << (filename + ext) << "..." << endl;
    ofstream out_file((filename + ext).c_str(), ios::binary);
    if (!out_file) {                                    // Check if file opened successfully
        cerr << "Error opening file " << (filename + ext) << " for writing!" << endl;
        exit(1);
    }
    // File header is P6 for binary pix map, x_size y_size, max pixel value per RGB component
    if (USE_COLOR) {
        out_file << "P6\n" << image_size << " " << image_size << "\n" << (MAX_PIXEL_LEVELS - 1) << "\n";
    } else {
        out_file << "P5\n" << image_size << " " << image_size << "\n" << (MAX_PIXEL_LEVELS - 1) << "\n";
    }

    // Compute histogram
    int totalCount = 0;
    for (int j = 0; j < image_size; ++j) {
        for (int i = 0; i < image_size; ++i) {
            if (image[j][i] > 0) {
                ++histogram[image[j][i]];
                ++totalCount;
            }
        }
    }
    // Compute cumulative histogram
    int count = 0;
    for (int i = 1; i < max_iteration; ++i) {
        count += histogram[i];
        histogram[i] = count;
    }
    // Normalize using histogram
    for (int j = 0; j < image_size; ++j) {
        for (int i = 0; i < image_size; ++i) {
            if (USE_COLOR) {
                RGB rgb = BLACK;
                if (image[j][i] != 0) {
#ifdef USE_RAINBOW
                    rgb = level_2_RGB(histogram[image[j][i]] / (double) (totalCount + 1));
#else
                    rgb = blend_colors(histogram[image[j][i]] / (double) (totalCount + 1), color1, color2);
#endif
                }
#if MAX_PIXEL_LEVELS > 256
                out_file.write(reinterpret_cast<const char*>(&rgb.R) + 1, 1);
                out_file.write(reinterpret_cast<const char*>(&rgb.R), 1);
                out_file.write(reinterpret_cast<const char*>(&rgb.G) + 1, 1);
                out_file.write(reinterpret_cast<const char*>(&rgb.G), 1);
                out_file.write(reinterpret_cast<const char*>(&rgb.B) + 1, 1);
                out_file.write(reinterpret_cast<const char*>(&rgb.B), 1);
#else
                out_file.write(reinterpret_cast<const char*>(&rgb.R), 1);
                out_file.write(reinterpret_cast<const char*>(&rgb.G), 1);
                out_file.write(reinterpret_cast<const char*>(&rgb.B), 1);
#endif
            } else {
#if MAX_PIXEL_LEVELS > 256
                uint16_t level = 0;
                if (image[j][i] != 0)
                    level = color1.R + histogram[image[j][i]] / (double) (totalCount + 1) * (color2.R - color1.R + 1);

                out_file.write(reinterpret_cast<const char*>(&level) + 1, 1);
                out_file.write(reinterpret_cast<const char*>(&level), 1);
#else
                uint8_t level = 0;
                if (image[j][i] != 0)
                    level = color1.R + histogram[image[j][i]] / (double) (totalCount + 1) * (color2.R - color1.R + 1);
                out_file.write(reinterpret_cast<const char*>(&level), 1);
#endif
            }
        }
    }
    out_file.close();
}

/**
* @brief Main function.
* 
* The entry point of the program.
*
* Syntax: mandelbrot <x> <y> <radius> [--max_iterations=n] [--image_size=n] 
*   [--filename=name] [--high_color=color] [--low_color=color] [--decimal_digits=n]
*
* x & y specify the coordinates of the Mandelbrot Set to display at the center of the image.
*
* radius is the radius in Mandelbrot Set coordinates out to the extreme edges of the image.
*
* max_iterations is the number of iterations maximum per pixel to determine if it is in the 
* set.  The max_iterations will need to get larger as radius gets smaller.
*
* image_size is the size in pixels for the width and height of the image.
*
* high_color and low_color are hexadecimal color values.
*
* filename is the PGM or PPM file that will be created with the image.
* 
* decimal_digits sets the accuracy of floating point computations when using GMP.
*
* @return int Returns 0 on successful execution.
*/
int main(int argc, char** argv) {
    CLI::App app{"Program to generate Mandelbrot images."};
    app.set_help_flag("", "");
    // use CLI11 to process command line options.
    vector<string> coords;
    app.add_option("coords", coords, "Location to plot (x y radius)")
        ->expected(3);

    // Optional arguments
    string filename = "image";
    int max_iteration = 10000;
    int image_size = 500;
    string high_color = "#FFF";
    string low_color = "#000";
    bool show_help = false;
    int decimal_digits = 16;

    app.add_flag("--help", show_help, "Show help");
    app.add_option("-f,--filename", filename, "Base filename (no extension)");
    app.add_option("-m,--max_iterations", max_iteration, "Max. iterations");
    app.add_option("-s,--image_size", image_size, "Image size in pixels");
    app.add_option("-h,--high_color", high_color, "High color hex code (#xxx or #xxxxxx)");
    app.add_option("-l,--low_color", low_color, "Low color hex code (#xxx or #xxxxxx)");
    app.add_option("-d,--decimal_digits", decimal_digits, "Decimal digits of accuracy to use");

    CLI11_PARSE(app, argc, argv);

    if (show_help) {
        cout << app.help() << endl;
        return 0;
    }

    // Get the three required coordinates to plot
    double_type x, y, radius;
    if (coords.size() == 3) {
#ifdef MPREAL
        DIGITS = decimal_digits;
        mpreal::set_default_prec(PRECISION);    // Set default precision in bits
        x = double_type(coords[0]);
        y = double_type(coords[1]);
        radius = double_type(coords[2]);
#elif defined __FLOAT128__
        x = strtoflt128(coords[0], nullptr);
        y = strtoflt128(coords[1], nullptr);
        radius = strtoflt128(coords[2], nullptr);
#elif defined LONG_DOUBLE
        x = stold(coords[0]);
        y = stold(coords[1]);
        radius = stold(coords[2]);
#else
        x = stod(coords[0]);
        y = stod(coords[1]);
        radius = stod(coords[2]);
#endif
    } else {
        cout << "Arguments x, y and radius required but not found!" << endl;
        return 1;
    }

    // Check the image size
    if (image_size > sqrt(numeric_limits<int>::max())) {
        cout << "Image size can not exceed " << sqrt(numeric_limits<int>::max()) << endl;
        return 1;
    }
    
    // Get the low and high colors
    RGB high_RGB = RGB(high_color);
    RGB low_RGB = RGB(low_color);
    if (high_RGB.isGray() && low_RGB.isGray()) USE_COLOR = false;

    // Compute limits of image
    double_type x_min = x - radius;
    double_type x_max = x + radius;
    double_type y_min = y - radius;
    double_type y_max = y + radius;

    // Output settings
    cout << "Size of floating point: " << sizeof(double_type) << endl;
    cout << "Floating point decimal digits: " << DOUBLE_DIGITS << endl;
    cout << "Floating point precision: " << PRECISION << " bits" << endl;
    cout << "Maximum Iterations: " << max_iteration << endl;
    cout << "High escape color: " << high_color << endl;
    cout << "Low escape color: " << low_color << endl;
    cout << (USE_COLOR ? "Using color" : "Using grayscale") << endl;
    cout << "Image Bounds: " << setprecision(DOUBLE_DIGITS) << x_min << ", " << x_max << ", " << y_min << ", " << y_max << endl;
    cout << "Image Size: " << image_size << " x " << image_size << endl;
    
    // Allocate an array for the row pointers.
    image = make_unique<unique_ptr<int32_t[]>[]>(image_size);
    if (image == nullptr) {                             // Check for failure
        cerr << "Memory allocation failed!" << endl;
        exit(1);
    }

    // Do a simple check to see if our doubles have enough precision
    double_type x0 = x_max - ((x_max - x_min) / image_size) * (1 + double_type(0.5));
    double_type x1 = x_max - ((x_max - x_min) / image_size) * (2 + double_type(0.5));
    if (x0 == x1) 
        cout << "Double type may not have enough precision." << endl;

    // Get the number of CPU cores.
    unsigned int num_cpus = thread::hardware_concurrency();
    if (num_cpus == 0) {
        cout << "Could not determine number of CPUs.  Using 4.\n";
        num_cpus = 4;
    } else {
        cout << "Number of CPU cores: " << num_cpus << endl;
    }

    // Create the threads
    vector<thread> threads;
    for (int i = 0; i < num_cpus; ++i) {
        threads.emplace_back(compute_mandelbrot, x_min, x_max, y_min, y_max, max_iteration, image_size);
    }
    // Wait for all threads to complete
    for (auto& t : threads) {
        t.join();  
    }

    // Generate the output file
    output_file(filename, max_iteration, image_size, low_RGB, high_RGB);
    return 0;
}