Simple Mandelbrot Set Program

/**
* @file mandelbrot2.cpp
* @brief Program to compute Mandelbrot Set images.
* 
* This program can compute images of the Mandelbrot Set in its
* entirety or close-ups of boundary areas.
*
* Requires:
* Linux<br>
* get cli11.hpp at https://github.com/CLIUtils/CLI11/releases
*
* macOS<br>
* brew install cli11<br>
*
* Build:
* Linux<br>
* put CLI.hpp in /usr/local/include <br>
*
* macOS<br>
* Add to .zshrc<br>
* export CPPFLAGS="-I/opt/homebrew/include"<br>
* export LDFLAGS="-L/opt/homebrew/lib"
*
* Compile with<br>
* Linux<br>
* g++ --std=gnu++20 -o mandelbrot2 mandelbrot2.cpp

* macOS<br>
* g++ $CPPFLAGS $LDFLAGS --std=gnu++20 -o mandelbrot2 mandelbrot2.cpp
*
* Execute:
* mandelbrot2 real imaginary radius [--filename=basename]
*     [--image_size=pixels] [--max_iterations=count]
*
* @author Richard Lesh
* @date 2025-03-12
* @version 1.0
*/

#include <CLI/CLI.hpp>
#include <complex>
#include <fstream>
#include <string>
#include <vector>

#define MAX_PIXEL_LEVELS 256            // 256 for 8-bit

using namespace std;

/**
* @brief Computes the Unit in the Last Place (ULP) of a double-precision value.
*
* This function calculates the distance between the given value @p x and the next
* representable double greater than @p x in the direction of positive infinity.
* The result gives a measure of the precision at the value @p x.
*
* @param x The input double value.
* @return The ULP of @p x (i.e., the difference to the next larger representable double).
*
* @note If @p x is NaN or infinity, the result may be undefined.
*/
inline double ULP(double x) {
    double next = nextafter(x, numeric_limits<double>::infinity());
    double ulp = next - x;
    return ulp;
}

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

/**
* @brief Computes the squared magnitude (norm) of a complex number without
* taking the square root.
*
* This function returns the sum of the squares of the real and imaginary parts
* of the complex number @p z.  It avoids the overhead of computing a square
* root, making it useful for performance-sensitive contexts like Mandelbrot
* set calculations where only relative magnitude comparisons are needed.
*
* @tparam T The numeric type (e.g., float, double, mpreal).
* @param z The complex number whose squared magnitude is to be computed.
* @return The squared magnitude of @p z, equivalent to |z|^2.
*
* @see std::norm may be used as a standard alternative, but may perform
* additional checks.
*/

template<typename T>
inline T norm_squared(complex<T> z) {
    return z.real() * z.real() + z.imag() * z.imag();
}

// Number of pixels that escape
long num_escaped = 0L;
// Number of cycles detected
long num_cycles = 0L;
// Number of pixels that reach max. iterations
long num_max_iter = 0L;

/**
* @brief Compute the number of iterations to escape.
*
* Function computes the number of iterations of the formula
* z(n+1) = z(n)^2 + c that can be executed before the value of |z|
* exceeds the escape limit radius of 2.0.
*
* The function uses Floyd’s Cycle Detection (Tortoise 
* and Hare) algorithm to detect cycles.

* @param c Complex point used to perform the calculation.
* @param max_iteration Number of iterations, if exceeded, when we will declare
*                      as not escaping and therefore in the Mandelbrot set.
* @return Value [0, max_iteration) representing the escape iterations.
*                      0 is considered as not escaping.
*/

int compute_escape(complex<double> c, int max_iteration) {
    complex<double> tortoise{0., 0.};
    complex<double> hare{0., 0.};

    int i = 0;
    while (i < max_iteration && norm_squared(tortoise) < 4.0) {
        tortoise = tortoise * tortoise + c;
        hare = hare * hare + c;
        hare = hare * hare + c; // hare moves twice

        double ulps = ULP(std::norm(tortoise));
        ulps = 9 * ulps * ulps;
        if (norm_squared(tortoise - hare) < ulps) {
            ++num_cycles;
            return 0L; // cycle detected
        }
        ++i;
    }

    if (i == max_iteration) {   // Return 0 (black) if we are in the set
        ++num_max_iter;
        return 0L;
    } else {
        ++num_escaped;
        return i;
    }
}

/**
* @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 square 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_per_pixel = (x_max - x_min) / image_size;
    double height_per_pixel = (y_max - y_min) / image_size;
    int modulus = ceil(image_size / 40.);     // Used to control progress output
    if (row >= 0 && row < image_size) {
        image[row] = make_unique<int32_t[]>(image_size);
        // Calculate the imaginary coordinate of the middle of the pixel row
        double y0 = y_max - height_per_pixel * (row + 0.5);
        for (int i = 0; i < image_size; ++i) {
            // Calculate the real coordinate of the middle of the pixel
            double x0 = width_per_pixel * (i + 0.5) + x_min;
            // Compute the escape number of iterations [0, max_iteration)
            int level = compute_escape(complex(x0, y0), max_iteration);
            // Store the escape iterations in the image array.
            image[row][i] = level;
        }
        // Output periodic progress
        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 used to compute each point.
* @param image_size Size, in pixels, of the image to save.
*/

void output_file(const string &filename, int max_iteration, int image_size) {
    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
    out_file << "P5\n" << image_size << " " << image_size << "\n" << (MAX_PIXEL_LEVELS - 1) << "\n";

    // Find the maximum iteration count
    int max = 0;
    for (int j = 0; j < image_size; ++j) {
        for (int i = 0; i < image_size; ++i) {
            if (image[j][i] > max) max = image[j][i];
        }
    }
    max_iteration = max + 1;
    // Compute bucket size for histogram
    int max_iteration_bits = static_cast<int>(ceil(log(max_iteration)/log(2.0)));
    int bucket_shift_bits = max_iteration_bits > 20 ? max_iteration_bits - 20 : 0;
    vector<int>::size_type histogram_size = (max_iteration >> bucket_shift_bits) + 1;
    // Compute histogram
    vector<int> histogram(histogram_size, 0);
    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[vector<int>::size_type(image[j][i] >> bucket_shift_bits)];
                ++totalCount;
            }
        }
    }
    // Compute cumulative histogram
    int count = 0;
    for (vector<int>::size_type i = 1; i < histogram_size; ++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 = floor(histogram[vector<int>::size_type(image[j][i] >> bucket_shift_bits)]
                / static_cast<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 real imaginary radius [--filename=basename]
*     [--image_size=n] [--max_iterations=n]
* 
* real & imaginary 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 maximum number of iterations per pixel allowed to
* determine if it is in the set.  The max_iterations will generally need to get
* larger as radius gets smaller.
*
* image_size is the size in pixels for the width and height of the square image.
*
* filename is the basename of the image file.  It defaults to "image".
*
* @return int Returns 0 on successful execution.
*/
int main(int argc, char** argv) {
    CLI::App app{"Program to generate Mandelbrot Set images."};
    app.set_help_flag("", "");
    argv = app.ensure_utf8(argv);
    // 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 = 512;
    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");

    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[vector<string>::size_type(0)]);
        y = stod(coords[vector<string>::size_type(1)]);
        radius = stod(coords[vector<string>::size_type(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 " << floor(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);

    // Do a simple check to see if our doubles have enough precision
    if ((x_max - x_min) / image_size < ULP(x) ||
        (y_max - y_min) / image_size < ULP(y))
        cout << "Doubles may not have enough precision!" << endl;

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

    cout << endl;
    cout << "Number of escaped pixels: " << num_escaped << endl;
    cout << "Number of cycles detected: " << num_cycles << endl;
    cout << "Number of max iteration pixels: " << num_max_iter << endl;

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