//
//  main.cpp
//  subpixel_edge
//
//  Created by Mehdi Stapleton on 2015-05-18.
//  Copyright (c) 2015 Mehdi Stapleton. All rights reserved.
//

#include "subpixel_edge.h"

//---------------------------------------------------------------
//                          DESCRIPTION
//---------------------------------------------------------------

/*Mex Function Inputs:
 1)X gradient image     Desc:X-component (column-wise) of gradient.
                             Should be preprocessed for removing any 'small' gradients
                        Type:NxM matrix
 2)Y gradient image     Desc:Y-component (row-wise) of gradient.
                             Should be preprocessed for removing any 'small' gradients
                        Type:NxM matrix
 */

/*Mex Function Outputs:
 1)Edge-magnitude image Desc:Magnitude of non-maximally-suppressed edges (|gradient|'s)
                        Type:NxM matrix
 2)Delta X image        Desc:Subpixel offset (column-wise) from integer edge pixel location
                        Type:NxM matrix
 3)Delta Y image        Desc:Subpixel offset (row-wise) from integer edge pixel location
                        Type:NxM matrix
 */

/*Background
 >performs non-maximum suppression (with quadratic interpolation on real gradient) and then subsequently applies
  quadratic interpolation to locate maximum gradient value in direction of gradient (computing subpixel offsets)
 >see: Devernay, Frédéric. 
 "A non-maxima suppression method for edge detection with sub-pixel accuracy." (1995). 
 */


#ifndef XCODE_COMPILATION

//-------------------------------------------------------------------------
//              MEX File Generator Function Entry-Point
//-------------------------------------------------------------------------

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    /* Macros for the ouput and input arguments */
#define XGRAD_IMAGE_IN      prhs[0]
#define YGRAD_IMAGE_IN      prhs[1]
    
#define EDGE_IMAGE_OUT     plhs[0]
#define XDELTA_IMAGE_OUT   plhs[1]
#define YDELTA_IMAGE_OUT   plhs[2]
    
    matrix_t *xGradImgIn;
    matrix_t *yGradImgIn;
    
    matrix_t *edgeImgOut;
    matrix_t *xDeltaImgOut;
    matrix_t *yDeltaImgOut;
    
    Image<grad_t> xGradImgOcl;
    Image<grad_t> yGradImgOcl;
    
    Image<grad_t>     edgeImgOcl;
    Image<subpixel_t> xDeltaImgOcl;
    Image<subpixel_t> yDeltaImgOcl;
    
    mwSize        nImgDims;
    const mwSize* imgDims;
    size_t        nRows,nCols;
    
    if(nrhs != 2) /* Check the number of input arguments */
    {
        LOG_ERROR("Wrong number of input arguments.");
    }
    else if(nlhs != 3) /* Check the number of output arguments */
    {
        LOG_ERROR("Too many output arguments.");
    }
    
    if(!IS_REAL_2D_FULL_DOUBLE(XGRAD_IMAGE_IN) || !IS_REAL_2D_FULL_DOUBLE(YGRAD_IMAGE_IN))
    {
        LOG_ERROR("Input (x/y)-gradient images must be a real 2D full double arrays.");
    }
    
    nImgDims = mxGetNumberOfDimensions(XGRAD_IMAGE_IN);
    imgDims  = mxGetDimensions(XGRAD_IMAGE_IN);
    nRows    = imgDims[0];
    nCols    = imgDims[1];
    
    if( nRows <= 0 || nCols <= 0)
    {
        LOG_ERROR("Invalid number of rows or columns");
    }
    
    LOG_DEBUG("Received color image with dimensions " << nRows << "x" << nCols);
    
    xGradImgIn = mxGetPr(XGRAD_IMAGE_IN);
    yGradImgIn = mxGetPr(YGRAD_IMAGE_IN);
    
    EDGE_IMAGE_OUT   = mxCreateNumericArray(nImgDims, imgDims, mxDOUBLE_CLASS, mxREAL);
    XDELTA_IMAGE_OUT = mxCreateNumericArray(nImgDims, imgDims, mxDOUBLE_CLASS, mxREAL);
    YDELTA_IMAGE_OUT = mxCreateNumericArray(nImgDims, imgDims, mxDOUBLE_CLASS, mxREAL);
    
    edgeImgOut   = (matrix_t*)mxGetData(EDGE_IMAGE_OUT);
    xDeltaImgOut = (matrix_t*)mxGetData(XDELTA_IMAGE_OUT);
    yDeltaImgOut = (matrix_t*)mxGetData(YDELTA_IMAGE_OUT);
    
    xGradImgOcl.width  = nCols;
    xGradImgOcl.height = nRows;
    yGradImgOcl.width  = nCols;
    yGradImgOcl.height = nRows;
    
    edgeImgOcl.width    = nCols;
    edgeImgOcl.height   = nRows;
    xDeltaImgOcl.width  = nCols;
    xDeltaImgOcl.height = nRows;
    yDeltaImgOcl.width  = nCols;
    yDeltaImgOcl.height = nRows;
    
    //(matlab column-major -> opencl row-major)
    for(index_t r=0; r<nRows; ++r)
    {
        for(index_t c=0; c<nCols; ++c)
        {
            xGradImgOcl.pixels.push_back((grad_t)xGradImgIn[r+nRows*c]);
            yGradImgOcl.pixels.push_back((grad_t)yGradImgIn[r+nRows*c]);
        }
    }
    
    //initialize to all-zeros output images
    edgeImgOcl.zero_fill();
    xDeltaImgOcl.zero_fill();
    yDeltaImgOcl.zero_fill();
    
    //perform GPU-accelerated operation
    subpixel_edge_gpu(xGradImgOcl, yGradImgOcl, edgeImgOcl, xDeltaImgOcl, yDeltaImgOcl);
                          
    //(opencl row-major -> matlab column-major)
    for(index_t r=0; r<nRows; ++r)
    {
        for(index_t c=0; c<nCols; ++c)
        {
            edgeImgOut[r+nRows*c]   = (matrix_t)edgeImgOcl.pixels[c+nCols*r];
            xDeltaImgOut[r+nRows*c] = (matrix_t)xDeltaImgOcl.pixels[c+nCols*r];
            yDeltaImgOut[r+nRows*c] = (matrix_t)yDeltaImgOcl.pixels[c+nCols*r];
        }
    }
    
    
}//end of mexFunction

#endif//XCODE_COMPILATION