//a example mex c-file for use as a basic template
//[C]=cMEXexample(A,B)


//mex -O cMEXexample.c

//basic includes, others may be needed depending on application
#include "mex.h"
#include "math.h"

// This function is mandatory.  It takes the place of 'main' in c.
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
    
    double *pA, *pB, *pC; // pointers to incoming and outgoing data
    int N; //number of pixles
    int *pSizeA, *pSizeB; // size array
    int dims[2]; // output size array
    int Xsize, Ysize;
    int Adims, Bdims;
    int xx, yy;
    
    /* Check for proper number of arguments. */
    if (nrhs != 2) {
        mexErrMsgTxt("two inputs required.");
    }
    
    /* Assign pointers to each input. */
    pA=mxGetPr(prhs[0]);
    pB=mxGetPr(prhs[1]);
    
    //get dimensions and sizes
    Adims=mxGetNumberOfDimensions(prhs[0]);
    Bdims=mxGetNumberOfDimensions(prhs[1]);
    
    pSizeA = mxGetDimensions(prhs[0]);
    pSizeB = mxGetDimensions(prhs[1]);
    
    // this writes to the matlab command line.
    mexPrintf("size A:  %d  X  %d\n", pSizeA[0], pSizeA[1]);
    mexPrintf("size B:  %d  X  %d\n", pSizeB[0], pSizeB[1]);
    
    
    //check dimensions
    if (Adims != 2) {
        mexErrMsgTxt("A must be 2D.");
    }
    
    if (Bdims != 2) {
        mexErrMsgTxt("B must be 2D.");
    }
    
    //check sizes
    if (pSizeA[0] != pSizeB[0]) {
        mexErrMsgTxt("A and B must be same size.");
    }
    
    if (pSizeA[1] != pSizeB[1]) {
        mexErrMsgTxt("A and B must be same size.");
    }
    
    //Note that we did not check data type and assumed it was 'double'.
    //This is an example of bad programming (on purpose this time)
    
    Xsize=pSizeA[0];
    Ysize=pSizeA[1];
    
    /* Create matrix for the return argument. */
    dims[0]=(int)pSizeA[0];
    dims[1]=(int)pSizeA[1];
    
    plhs[0] = mxCreateDoubleMatrix((int)pSizeA[0], ((int)pSizeA[1]), mxREAL);
    
    pC = mxGetPr(plhs[0]);
     
    //write into memory assigned to output matrix
    for (xx=0;xx<Xsize;xx++)
        for (yy=0;yy<Ysize;yy++) {
        pC[yy*Xsize+xx]=pA[yy*Xsize+xx]*pB[yy*Xsize+xx];
        
        }
    
    
    
}