DCTPar.cpp

#include <iostream>
#include <vector>
#include <string>
#include <cmath>
#include <algorithm>
#include <chrono>
#include <opencv2/opencv.hpp>
#include <omp.h>

using namespace std;
using namespace cv;

static const float PI = 3.14159265358979323846f;


///////////////////////////////////
//// DCT TABLES STRUCTURE /////////
///////////////////////////////////
struct DCTTables {
    int N;
    vector<float> alpha;
    vector<float> cosTable;
    DCTTables(int N_) : N(N_), alpha(N_), cosTable(N_ * N_) {
        for (int k = 0; k < N; ++k) {
            alpha[k] = (k == 0) ? sqrtf(1.0f / N) : sqrtf(2.0f / N);
        }
        for (int k = 0; k < N; ++k) {
            for (int n = 0; n < N; ++n) {
                cosTable[k * N + n] = cosf(PI * (2.0f * n + 1.0f) * k / (2.0f * N));
            }
        }
    }
};


////////////////////////////////
//// MATRIX OPERATIONS /////////
////////////////////////////////
vector<vector<float>> subtractMatrices(const vector<vector<float>>& a, const vector<vector<float>>& b) {
    int H = a.size();
    int W = a[0].size();

    vector<vector<float>> result(H, vector<float>(W));

    for (int r = 0; r < H; r++) {
        for (int c = 0; c < W; c++) {
            result[r][c] = a[r][c] - b[r][c];
        }
    }

    return result;
}


///////////////////////////////
/////// IMAGE LOADING /////////
///////////////////////////////
vector<vector<float>> loadImageFromDisk(const string& filename) {

    //// GREYSCALE ////
    Mat img = imread(filename, IMREAD_GRAYSCALE);

    if (img.empty()) {
        cerr << "Errore: impossibile caricare l'immagine " << filename << endl;
        return {};
    }

    int H = img.rows;
    int W = img.cols;

    //// CONVERT TO VECTOR<VECTOR<FLOAT>> ////
    vector<vector<float>> result(H, vector<float>(W));

    for (int r = 0; r < H; r++) {
        for (int c = 0; c < W; c++) {
            result[r][c] = (float)img.at<uchar>(r, c);
        }
    }

    return result;
}

///////////////////////////////////
/////// BLOCK EXTRACTION //////////
///////////////////////////////////
vector<vector<float>> extractBlock(const vector<vector<float>>& img, int startRow, int startCol, int blockSize) {

    vector<vector<float>> block(blockSize, vector<float>(blockSize, 0.0f));
    int imgHeight = img.size();
    int imgWidth = img[0].size();

    //// COPY BLOCK DATA ////
    for (int r = 0; r < blockSize; r++) {
        for (int c = 0; c < blockSize; c++) {
            if (startRow + r < imgHeight && startCol + c < imgWidth) {
                block[r][c] = img[startRow + r][startCol + c];
            }
        }
    }

    return block;
}

///////////////////////////////
///// 2D DCT WITH TABLES //////
///////////////////////////////
vector<vector<float>> dct2d(const vector<vector<float>>& img, const DCTTables& tables) {
    int H = img.size();
    int W = img[0].size();

    vector<vector<float>> out(H, vector<float>(W));

    //// PRECALCULATED TABLES ////
    const float* alphaArr = tables.alpha.data();
    const float* cosArr = tables.cosTable.data();

    //// DCT COMPUTATION ////
    for (int u = 0; u < H; u++) {
        for (int v = 0; v < W; v++) {
            float sum = 0.0f;

            for (int x = 0; x < H; x++) {
                float cos_u = cosArr[u * tables.N + x];
                for (int y = 0; y < W; y++) {
                    float cos_v = cosArr[v * tables.N + y];
                    sum += img[x][y] * cos_u * cos_v;
                }
            }

            //// STORE RESULT ////
            out[u][v] = alphaArr[u] * alphaArr[v] * sum;
        }
    }

    return out;
}

//////////////////////////////
//// 2D IDCT WITH TABLES ////
//////////////////////////////
vector<vector<float>> idct2d(const vector<vector<float>>& X, const DCTTables& tables) {
    int H = X.size();
    int W = X[0].size();

    vector<vector<float>> out(H, vector<float>(W));

    //// PRECALCULATED TABLES ////
    const float* alphaArr = tables.alpha.data();
    const float* cosArr = tables.cosTable.data();

    //// IDCT COMPUTATION ////
    for (int x = 0; x < H; x++) {
        for (int y = 0; y < W; y++) {
            float sum = 0.0f;

            for (int u = 0; u < H; u++) {
                float alpha_u = alphaArr[u];
                float cos_u = cosArr[u * tables.N + x];
                for (int v = 0; v < W; v++) {
                    float alpha_v = alphaArr[v];
                    float cos_v = cosArr[v * tables.N + y];
                    sum += alpha_u * alpha_v * X[u][v] * cos_u * cos_v;
                }
            }

            //// STORE RESULT ////
            out[x][y] = sum;
        }
    }

    return out;
}

//////////////////////////////
//// PROCESS IMAGE BLOCKS ////
//////////////////////////////
void processImageBlocks(const vector<vector<float>>& img,
                       vector<vector<float>>& dctComplete,
                       vector<vector<float>>& reconstructed,
                       int blockSize,
                       const DCTTables& tables) {

    int imgHeight = img.size();
    int imgWidth = img[0].size();

    int totalBlocks = ((imgHeight + blockSize - 1) / blockSize) *
                     ((imgWidth + blockSize - 1) / blockSize);
    int blocksProcessed = 0;

#pragma omp parallel default(none) shared(img, dctComplete, reconstructed, imgHeight, imgWidth, blockSize, blocksProcessed, tables)
    {
        //// LOCAL THREAD BUFFER ////
        vector<vector<float>> block(blockSize, vector<float>(blockSize));

#pragma omp for collapse(2) schedule(dynamic) nowait
        for (int startRow = 0; startRow < imgHeight; startRow += blockSize) {
            for (int startCol = 0; startCol < imgWidth; startCol += blockSize) {

                int actualH = min(blockSize, imgHeight - startRow);
                int actualW = min(blockSize, imgWidth - startCol);

                //// COPY DATA TO LOCAL BUFFER ////
                for (int r = 0; r < actualH; r++) {
                    for (int c = 0; c < actualW; c++) {
                        block[r][c] = img[startRow + r][startCol + c];
                    }
                }

                //// PADDING WITH ZEROS ////
                for (int r = actualH; r < blockSize; r++) {
                    for (int c = 0; c < blockSize; c++) {
                        block[r][c] = 0.0f;
                    }
                }
                for (int r = 0; r < actualH; r++) {
                    for (int c = actualW; c < blockSize; c++) {
                        block[r][c] = 0.0f;
                    }
                }

                //// APPLY DCT 2D ////
                auto X = dct2d(block, tables);

                //// APPLY IDCT 2D ////
                auto rec = idct2d(X, tables);

                //// COPY RESULTS TO COMPLETE MATRICES ////
                for (int r = 0; r < actualH; r++) {
                    for (int c = 0; c < actualW; c++) {
                        dctComplete[startRow + r][startCol + c] = X[r][c];
                        reconstructed[startRow + r][startCol + c] = rec[r][c];
                    }
                }

                #pragma omp atomic
                blocksProcessed++;
            }
        }
    }

}

///////////////////////
//// MAIN FUNCTION ////
///////////////////////
int main(int argc, char** argv) {

    //// THREAD CONFIGURATION ////
    int numThreads = omp_get_max_threads();
    if (argc >= 3) {
        numThreads = atoi(argv[2]);
        omp_set_num_threads(numThreads);
    }


    //// IMAGE LOAD ////
    string imagePath;
    if (argc >= 2) {
        imagePath = argv[1];
    } else {
        imagePath = "/home/lapemaya/PycharmProjects/SitoBorsa/lora/flux-controlnet-canny-Lora1.png";

    }
    auto img = loadImageFromDisk(imagePath);
    if (img.empty()) {
        cerr << "Errore: impossibile caricare l'immagine\n";
        return 1;
    }

    //// IMAGE VARIABLES ////
    int imgHeight = img.size();
    int imgWidth = img[0].size();
    int blockSize = 128;

    int numBlocksH = (imgHeight + blockSize - 1) / blockSize;
    int numBlocksW = (imgWidth + blockSize - 1) / blockSize;

    //// MEMORY ALLOC ////
    vector<vector<float>> dctComplete(imgHeight, vector<float>(imgWidth, 0.0f));
    vector<vector<float>> reconstructed(imgHeight, vector<float>(imgWidth, 0.0f));

    //// DCT TABLES ////
    DCTTables tables(blockSize);

    //// TIMER START ////
    auto start = chrono::high_resolution_clock::now();

    //// PROCESS IMAGE BLOCKS ////
    processImageBlocks(img, dctComplete, reconstructed, blockSize, tables);

    //// TIMER END ////
    auto end = chrono::high_resolution_clock::now();
    auto duration = chrono::duration_cast<chrono::milliseconds>(end - start);

    //// OUTPUT METRICS ////
    cout << "\n=== TEMPO DI ELABORAZIONE: " << duration.count() << " ms ===\n";
    cout << "Throughput: " << (numBlocksH * numBlocksW * 1000.0 / duration.count())
         << " blocchi/secondo\n\n";

    //// COMPUTE ERROR METRICS ////
    float maxError = 0.0f;
    float sumError = 0.0f;
    long long totalPixels = 0;

    for (int r = 0; r < imgHeight; r++) {
        for (int c = 0; c < imgWidth; c++) {
            float err = fabsf(reconstructed[r][c] - img[r][c]);
            maxError = max(maxError, err);
            sumError += err;
            totalPixels++;
        }
    }

    //// PRINT METRICS ////
    float avgError = sumError / totalPixels;
    cout << "=== METRICHE DI QUALITÀ ===\n";
    cout << "Errore medio di ricostruzione: " << avgError << "\n";
    cout << "Errore massimo di ricostruzione: " << maxError << "\n";
    cout << "Pixel totali processati: " << totalPixels << "\n\n";


    return 0;
}