SPMD.cpp

#include <iostream>
#include <vector>
#include <fstream>
#include <chrono>
#include <omp.h>
#include <cstdlib>

using namespace std;

// ==========================================================
// 8K RESOLUTION
// ==========================================================

const int WIDTH  = 7680;
const int HEIGHT = 4320;

const int MAX_ITER = 1000;

// ==========================================================
// PIXEL STRUCT
// ==========================================================

struct Pixel {
    unsigned char r, g, b;
};

// ==========================================================
// SAVE RGB IMAGE
// ==========================================================

void savePPM(const string& filename,
             const vector<Pixel>& image) {

    ofstream file(filename, ios::binary);

    file << "P6\n"
         << WIDTH << " "
         << HEIGHT << "\n255\n";

    for (const auto& pixel : image) {
        file.write(
            reinterpret_cast<const char*>(&pixel),
            3
        );
    }

    file.close();
}

// ==========================================================
// SAVE GRAYSCALE IMAGE
// ==========================================================

void saveGrayPPM(const string& filename,
                 const unsigned char* image,
                 size_t size) {

    ofstream file(filename, ios::binary);

    file << "P5\n"
         << WIDTH << " "
         << HEIGHT << "\n255\n";

    file.write(
        reinterpret_cast<const char*>(image),
        size
    );

    file.close();
}

// ==========================================================
// TASK A - MANDELBROT
// ==========================================================

void generateMandelbrot(vector<Pixel>& image) {

    double xmin = -2.5;
    double xmax = 1.0;
    double ymin = -1.0;
    double ymax = 1.0;

    #pragma omp parallel for schedule(guided)
    for (int y = 0; y < HEIGHT; y++) {

        for (int x = 0; x < WIDTH; x++) {

            double cr =
                xmin +
                (x / (double)WIDTH) *
                (xmax - xmin);

            double ci =
                ymin +
                (y / (double)HEIGHT) *
                (ymax - ymin);

            double zr = 0.0;
            double zi = 0.0;

            int iter = 0;

            while (zr * zr + zi * zi <= 4.0 &&
                   iter < MAX_ITER) {

                double temp =
                    zr * zr - zi * zi + cr;

                zi = 2.0 * zr * zi + ci;
                zr = temp;

                iter++;
            }

            Pixel pixel;

            if (iter == MAX_ITER) {

                pixel = {0, 0, 0};

            } else {

                unsigned char color =
                    (unsigned char)
                    (255.0 * iter / MAX_ITER);

                pixel.r = color;
                pixel.g = 255 - color;
                pixel.b = (color * 2) % 255;
            }

            image[y * WIDTH + x] = pixel;
        }
    }
}

// ==========================================================
// RGB -> GRAYSCALE
// ==========================================================

void convertToGray(
    const vector<Pixel>& input,
    unsigned char* gray)
{
    #pragma omp parallel for schedule(static)
    for (long long i = 0;
         i < (long long)(WIDTH * HEIGHT);
         i++) {

        gray[i] =
            (input[i].r +
             input[i].g +
             input[i].b) / 3;
    }
}

// ==========================================================
// TASK B - BLUR ROW (loop SIMD aislado)
// ==========================================================

__attribute__((noinline))
static void blurRow(
    const unsigned char* __restrict__ row_m2,
    const unsigned char* __restrict__ row_m1,
    const unsigned char* __restrict__ row_0,
    const unsigned char* __restrict__ row_p1,
    const unsigned char* __restrict__ row_p2,
    unsigned char* __restrict__ out_row)
{
    #pragma omp simd safelen(16) simdlen(16)
    for (int x = 2; x < WIDTH - 2; x++) {

        int sum =
            row_m2[x-2]*1 + row_m2[x-1]*4 + row_m2[x]*6 + row_m2[x+1]*4 + row_m2[x+2]*1 +
            row_m1[x-2]*4 + row_m1[x-1]*16+ row_m1[x]*24+ row_m1[x+1]*16+ row_m1[x+2]*4 +
            row_0 [x-2]*6 + row_0 [x-1]*24+ row_0 [x]*36+ row_0 [x+1]*24+ row_0 [x+2]*6 +
            row_p1[x-2]*4 + row_p1[x-1]*16+ row_p1[x]*24+ row_p1[x+1]*16+ row_p1[x+2]*4 +
            row_p2[x-2]*1 + row_p2[x-1]*4 + row_p2[x]*6 + row_p2[x+1]*4 + row_p2[x+2]*1;

        out_row[x] = sum / 256;
    }
}

// ==========================================================
// TASK B - GAUSSIAN BLUR
// ==========================================================

void gaussianBlurSIMD(const unsigned char* __restrict__ input,
                      unsigned char* __restrict__ output)
{
    const unsigned char* in =
        (const unsigned char*)__builtin_assume_aligned(input, 32);
    unsigned char* out =
        (unsigned char*)__builtin_assume_aligned(output, 32);

    #pragma omp parallel for schedule(static)
    for (int y = 2; y < HEIGHT - 2; y++) {

        blurRow(
            in + (y-2)*WIDTH,
            in + (y-1)*WIDTH,
            in + (y  )*WIDTH,
            in + (y+1)*WIDTH,
            in + (y+2)*WIDTH,
            out + y*WIDTH
        );
    }
}

// ==========================================================
// MAIN
// ==========================================================

int main() {

    omp_set_num_threads(2);

    cout << "=============================\n";
    cout << "OPENMP SIMD PROJECT\n";
    cout << "=============================\n";
    cout << "Threads: " << omp_get_max_threads() << "\n";

    // ======================================================
    // ARRAYS
    // ======================================================

    vector<Pixel> mandelbrot(WIDTH * HEIGHT);

    unsigned char* gray =
        (unsigned char*)aligned_alloc(32, WIDTH * HEIGHT);

    unsigned char* blurred =
        (unsigned char*)aligned_alloc(32, WIDTH * HEIGHT);

    // ======================================================
    // TASK A
    // ======================================================

    cout << "\nGenerating 8K Mandelbrot...\n";

    auto startA = chrono::high_resolution_clock::now();

    generateMandelbrot(mandelbrot);

    auto endA = chrono::high_resolution_clock::now();

    chrono::duration<double> timeA = endA - startA;

    cout << "Task A Time: " << timeA.count() << " seconds\n";

    savePPM("mandelbrot.ppm", mandelbrot);

    // ======================================================
    // RGB -> GRAY
    // ======================================================

    cout << "\nConverting to grayscale...\n";

    convertToGray(mandelbrot, gray);

    saveGrayPPM("gray.ppm", gray, WIDTH * HEIGHT);

    // ======================================================
    // TASK B
    // ======================================================

    cout << "\nApplying SIMD Gaussian Blur...\n";

    auto startB = chrono::high_resolution_clock::now();

    gaussianBlurSIMD(gray, blurred);

    auto endB = chrono::high_resolution_clock::now();

    chrono::duration<double> timeB = endB - startB;

    cout << "Task B Time: " << timeB.count() << " seconds\n";

    saveGrayPPM("blurred_simd.ppm", blurred, WIDTH * HEIGHT);

    cout << "\nDone.\n";

    free(gray);
    free(blurred);

    return 0;
}