idamnjanovic@60: /**************************************************************************
idamnjanovic@60:  *
idamnjanovic@60:  * File name: col2imstep.c
idamnjanovic@60:  *
idamnjanovic@60:  * Ron Rubinstein
idamnjanovic@60:  * Computer Science Department
idamnjanovic@60:  * Technion, Haifa 32000 Israel
idamnjanovic@60:  * ronrubin@cs
idamnjanovic@60:  *
idamnjanovic@60:  * Last Updated: 31.8.2009
idamnjanovic@60:  *
idamnjanovic@60:  *************************************************************************/
idamnjanovic@60: 
idamnjanovic@60: 
idamnjanovic@60: #include "mex.h"
idamnjanovic@60: 
idamnjanovic@60: 
idamnjanovic@60: /* Input Arguments */
idamnjanovic@60: 
idamnjanovic@60: #define	B_IN	 prhs[0]
idamnjanovic@60: #define N_IN   prhs[1]
idamnjanovic@60: #define SZ_IN  prhs[2]
idamnjanovic@60: #define S_IN   prhs[3]
idamnjanovic@60: 
idamnjanovic@60: 
idamnjanovic@60: /* Output Arguments */
idamnjanovic@60: 
idamnjanovic@60: #define	X_OUT	plhs[0]
idamnjanovic@60: 
idamnjanovic@60: 
idamnjanovic@60: void mexFunction(int nlhs, mxArray *plhs[], 
idamnjanovic@60: 		             int nrhs, const mxArray*prhs[])
idamnjanovic@60:      
idamnjanovic@60: { 
idamnjanovic@60:     double *x, *b, *s;
idamnjanovic@60:     mwSize sz[3], stepsize[3], n[3], ndims;
idamnjanovic@60:     mwIndex i, j, k, l, m, t, blocknum;
idamnjanovic@60:     
idamnjanovic@60:     
idamnjanovic@60:     /* Check for proper number of arguments */
idamnjanovic@60:     
idamnjanovic@60:     if (nrhs < 3 || nrhs > 4) {
idamnjanovic@60:       mexErrMsgTxt("Invalid number of input arguments."); 
idamnjanovic@60:     } else if (nlhs > 1) {
idamnjanovic@60:       mexErrMsgTxt("Too many output arguments."); 
idamnjanovic@60:     } 
idamnjanovic@60:     
idamnjanovic@60:     
idamnjanovic@60:     /* Check the the input dimensions */ 
idamnjanovic@60:     
idamnjanovic@60:     if (!mxIsDouble(B_IN) || mxIsComplex(B_IN) || mxGetNumberOfDimensions(B_IN)>2) {
idamnjanovic@60:       mexErrMsgTxt("B should be a double matrix.");
idamnjanovic@60:     }
idamnjanovic@60:     if (!mxIsDouble(N_IN) || mxIsComplex(N_IN) || mxGetNumberOfDimensions(N_IN)>2) {
idamnjanovic@60:       mexErrMsgTxt("Invalid output matrix size.");
idamnjanovic@60:     }
idamnjanovic@60:     ndims = mxGetM(N_IN)*mxGetN(N_IN);
idamnjanovic@60:     if (ndims<2 || ndims>3) {
idamnjanovic@60:       mexErrMsgTxt("Output matrix can only be 2-D or 3-D.");
idamnjanovic@60:     }
idamnjanovic@60:     if (!mxIsDouble(SZ_IN) || mxIsComplex(SZ_IN) || mxGetNumberOfDimensions(SZ_IN)>2 || mxGetM(SZ_IN)*mxGetN(SZ_IN)!=ndims) {
idamnjanovic@60:       mexErrMsgTxt("Invalid block size.");
idamnjanovic@60:     }
idamnjanovic@60:     if (nrhs == 4) {
idamnjanovic@60:       if (!mxIsDouble(S_IN) || mxIsComplex(S_IN) || mxGetNumberOfDimensions(S_IN)>2 || mxGetM(S_IN)*mxGetN(S_IN)!=ndims) {
idamnjanovic@60:         mexErrMsgTxt("Invalid step size.");
idamnjanovic@60:       }
idamnjanovic@60:     }
idamnjanovic@60:     
idamnjanovic@60:     
idamnjanovic@60:     /* Get parameters */
idamnjanovic@60:     
idamnjanovic@60:     s = mxGetPr(N_IN);
idamnjanovic@60:     if (s[0]<1 || s[1]<1 || (ndims==3 && s[2]<1)) {
idamnjanovic@60:       mexErrMsgTxt("Invalid output matrix size.");
idamnjanovic@60:     }
idamnjanovic@60:     n[0] = (mwSize)(s[0] + 0.01);
idamnjanovic@60:     n[1] = (mwSize)(s[1] + 0.01);
idamnjanovic@60:     n[2] = ndims==3 ? (mwSize)(s[2] + 0.01) : 1;
idamnjanovic@60:     
idamnjanovic@60:     s = mxGetPr(SZ_IN);
idamnjanovic@60:     if (s[0]<1 || s[1]<1 || (ndims==3 && s[2]<1)) {
idamnjanovic@60:       mexErrMsgTxt("Invalid block size.");
idamnjanovic@60:     }
idamnjanovic@60:     sz[0] = (mwSize)(s[0] + 0.01);
idamnjanovic@60:     sz[1] = (mwSize)(s[1] + 0.01);
idamnjanovic@60:     sz[2] = ndims==3 ? (mwSize)(s[2] + 0.01) : 1;
idamnjanovic@60:     
idamnjanovic@60:     if (nrhs == 4) {
idamnjanovic@60:       s = mxGetPr(S_IN);
idamnjanovic@60:       if (s[0]<1 || s[1]<1 || (ndims==3 && s[2]<1)) {
idamnjanovic@60:         mexErrMsgTxt("Invalid step size.");
idamnjanovic@60:       }
idamnjanovic@60:       stepsize[0] = (mwSize)(s[0] + 0.01);
idamnjanovic@60:       stepsize[1] = (mwSize)(s[1] + 0.01);
idamnjanovic@60:       stepsize[2] = ndims==3 ? (mwSize)(s[2] + 0.01) : 1;
idamnjanovic@60:     }
idamnjanovic@60:     else {
idamnjanovic@60:       stepsize[0] = stepsize[1] = stepsize[2] = 1;
idamnjanovic@60:     }
idamnjanovic@60:     
idamnjanovic@60:     if (n[0]<sz[0] || n[1]<sz[1] || (ndims==3 && n[2]<sz[2])) {
idamnjanovic@60:       mexErrMsgTxt("Block size too large.");
idamnjanovic@60:     }
idamnjanovic@60: 
idamnjanovic@60:     if (mxGetN(B_IN) != ((n[0]-sz[0])/stepsize[0]+1)*((n[1]-sz[1])/stepsize[1]+1)*((n[2]-sz[2])/stepsize[2]+1)) {
idamnjanovic@60:       mexErrMsgTxt("Invalid number of columns in B. Please use IM2COLSTEP to compute B.");
idamnjanovic@60:     }
idamnjanovic@60:     
idamnjanovic@60:     
idamnjanovic@60:     /* Create a matrix for the return argument */
idamnjanovic@60:     
idamnjanovic@60:     X_OUT = mxCreateNumericArray(ndims, n, mxDOUBLE_CLASS, mxREAL);
idamnjanovic@60:     
idamnjanovic@60:     
idamnjanovic@60:     /* Assign pointers */
idamnjanovic@60:     
idamnjanovic@60:     b = mxGetPr(B_IN);
idamnjanovic@60:     x = mxGetPr(X_OUT);
idamnjanovic@60:             
idamnjanovic@60:     
idamnjanovic@60:     /* Do the actual computation */
idamnjanovic@60:     
idamnjanovic@60:     blocknum = 0;
idamnjanovic@60:     
idamnjanovic@60:     /* iterate over all blocks */
idamnjanovic@60:     for (k=0; k<=n[2]-sz[2]; k+=stepsize[2]) {
idamnjanovic@60:       for (j=0; j<=n[1]-sz[1]; j+=stepsize[1]) {
idamnjanovic@60:         for (i=0; i<=n[0]-sz[0]; i+=stepsize[0]) {
idamnjanovic@60:           
idamnjanovic@60:           /* add single block */
idamnjanovic@60:           for (m=0; m<sz[2]; m++) {
idamnjanovic@60:             for (l=0; l<sz[1]; l++) {
idamnjanovic@60:               for (t=0; t<sz[0]; t++) {
idamnjanovic@60:                 (x+(k+m)*n[0]*n[1]+(j+l)*n[0]+i)[t] += (b + blocknum*sz[0]*sz[1]*sz[2] + m*sz[0]*sz[1] + l*sz[0])[t];
idamnjanovic@60:               }
idamnjanovic@60:             }
idamnjanovic@60:           }
idamnjanovic@60:           blocknum++;
idamnjanovic@60:           
idamnjanovic@60:         }
idamnjanovic@60:       }
idamnjanovic@60:     }
idamnjanovic@60:     
idamnjanovic@60:     return;
idamnjanovic@60: }