/**************************************************************************
 * Piotr's Image&Video Toolbox      Version 2.2
 * Copyright 2009 Piotr Dollar.  [pdollar-at-caltech.edu]
 * Please email me if you find bugs, or have suggestions or questions!
 * Licensed under the Lesser GPL [see external/lgpl.txt]
 *************************************************************************/
#include "mex.h"

#define arraysum(A, m, s, e, i) m=0; for(i=s; i<=e; i++) { m+=A[i]; };
int	maxi( int x, int y ) { return (x > y) ? x : y; };
int	mini( int x, int y ) { return (x < y) ? x : y; };

/**************************************************************************
 * Row non-linear seperable filter - sum (see nlfiltersep.m)
 *  x = nlfiltersep_sum( [ 1 9; 5 9; 0 0; 4 8; 7 3; 2 6], 1, 1 )
 *  y = [6 18; 6 18; 9 17; 11 11; 13 17; 9 9]; x-y
 * B(i,j) is the sum of A(i-r1:i+r2,j). It has the same dims as A.
 * This can be implemented effiicently because:
 *  B[i] = B[i-1] + A[i+r2] - A[i-r1-1];
 * Does not work for initial (r1+1) and final r2 values in each row.
 *************************************************************************/
void nlfiltersep_sum( const double *A, double *B, int r1, int r2, int mRows, int nCols ) {
  int i, row0, e, s, r, c; double m;
  for( c=0; c<nCols; c++ ) {
    row0 = mRows * c;
    /* leading border calculations */
    for(r=0; r<=mini(r1, mRows-1); r++) {
      e = mini( r+r2, mRows-1 );
      arraysum( A, m, row0, row0+e, i );
      B[r+row0] = m;
    }
    /* main caclulations */
    for(r=r1+1; r<mRows-r2; r++) {
      B[r+row0] = B[r+row0-1] + A[r+row0+r2] - A[r+row0-r1-1];
    }
    /* end border calculations */
    for(r=maxi(mRows-r2, 0); r<mRows; r++) {
      s = maxi( r-r1, 0 );
      arraysum( A, m, s+row0, mRows-1+row0, i );
      B[r+row0] = m;
    }
  }
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
  int mRows, nCols, r1, r2; double *A, *B;
  
  /* Error checking on arguments */
  if( nrhs!=3 ) mexErrMsgTxt("Three input arguments required.");
  if( nlhs>1 ) mexErrMsgTxt("Too many output arguments.");
  mRows = mxGetM(prhs[0]); nCols = mxGetN(prhs[0]);
  if(!mxIsDouble(prhs[0])) mexErrMsgTxt("Input array must be of type double.");
  
  /* extract inputs */
  A = (double*) mxGetData(prhs[0]);
  r1 = (int) mxGetScalar(prhs[1]);
  r2 = (int) mxGetScalar(prhs[2]);
  
  /* create outputs */
  plhs[0] = mxCreateDoubleMatrix(mRows, nCols, mxREAL );
  B = (double*) mxGetData(plhs[0]);
  
  /* Apply filter */
  nlfiltersep_sum( A, B, r1, r2, mRows, nCols );
}
