operator_idct_2d.cpp

// Software: Inverse Discrete Cosine Transform Operator 2D for JPEG
// Author: Hy Truong Son
// Position: PhD Student
// Institution: Department of Computer Science, The University of Chicago
// Email: sonpascal93@gmail.com, hytruongson@uchicago.edu
// Website: http://people.inf.elte.hu/hytruongson/
// Copyright 2016 (c) Hy Truong Son. All rights reserved.

#include <iostream>
#include <fstream>
#include <sstream>
#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <vector>
#include <set>
#include <iterator>
#include <algorithm>
#include <ctime>

#include "mex.h"

using namespace std;

void vector2matrix(double *input, int nRows, int nCols, double **output) {
    for (int i = 0; i < nRows; ++i) {
        for (int j = 0; j < nCols; ++j) {
            output[i][j] = input[j * nRows + i];
        }
    }
}

void matrix2vector(double **input, int nRows, int nCols, double *output) {
    for (int i = 0; i < nRows; ++i) {
        for (int j = 0; j < nCols; ++j) {
            output[j * nRows + i] = input[i][j];
        }
    }
}

double alpha(int u) {
    if (u == 0) {
        return 1.0 / sqrt(2.0);
    }
    return 1.0;
}

void mexFunction(int nOutputs, mxArray *output_pointers[], int nInputs, const mxArray *input_pointers[]) {
    if (nInputs != 2) {
        std::cerr << "The number of input parameters must be exactly 2 (Number of rows, and Number of columns)!" << std::endl;
        return;
    }

    // Number of rows - First input parameter
    int M = mxGetScalar(input_pointers[0]);
    
    // Number of columns - Second input parameter
    int N = mxGetScalar(input_pointers[1]);
    
    // Memory allocation
    /*
    double **block = new double* [M];
    for (int i = 0; i < M; ++i) {
        block[i] = new double [N];
    }
    */
    
    double **result = new double* [M * M];
    for (int i = 0; i < M * M; ++i) {
        result[i] = new double [N * N];
    }
    
    // Computation
    double c = 0.25;
    for (int x = 0; x < M; ++x) {
        for (int y = 0; y < N; ++y) {
            for (int u = 0; u < M; ++u) {
                for (int v = 0; v < N; ++v) {
                    double angle1 = (2.0 * x + 1.0) * (double)(u) * M_PI / (2.0 * M);
                    double angle2 = (2.0 * y + 1.0) * (double)(v) * M_PI / (2.0 * N);
                    result[x * M + u][y * N + v] = c * alpha(u) * alpha(v) * cos(angle1) * cos(angle2);
                }
            }
        }
    }
    
    // Return the result
    output_pointers[0] = mxCreateDoubleMatrix(M * M, N * N, mxREAL);
    matrix2vector(result, M * M, N * N, mxGetPr(output_pointers[0]));
}