Mercurial > hg > smallbox

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Problems/private/col2imstep.c	Thu Mar 25 14:05:27 2010 +0000
@@ -0,0 +1,146 @@
+/**************************************************************************
+ *
+ * File name: col2imstep.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 31.8.2009
+ *
+ *************************************************************************/
+
+
+#include "mex.h"
+
+
+/* Input Arguments */
+
+#define	B_IN	 prhs[0]
+#define N_IN   prhs[1]
+#define SZ_IN  prhs[2]
+#define S_IN   prhs[3]
+
+
+/* Output Arguments */
+
+#define	X_OUT	plhs[0]
+
+
+void mexFunction(int nlhs, mxArray *plhs[],
+		             int nrhs, const mxArray*prhs[])
+
+{
+    double *x, *b, *s;
+    mwSize sz[3], stepsize[3], n[3], ndims;
+    mwIndex i, j, k, l, m, t, blocknum;
+
+
+    /* Check for proper number of arguments */
+
+    if (nrhs < 3 || nrhs > 4) {
+      mexErrMsgTxt("Invalid number of input arguments.");
+    } else if (nlhs > 1) {
+      mexErrMsgTxt("Too many output arguments.");
+    }
+
+
+    /* Check the the input dimensions */
+
+    if (!mxIsDouble(B_IN) || mxIsComplex(B_IN) || mxGetNumberOfDimensions(B_IN)>2) {
+      mexErrMsgTxt("B should be a double matrix.");
+    }
+    if (!mxIsDouble(N_IN) || mxIsComplex(N_IN) || mxGetNumberOfDimensions(N_IN)>2) {
+      mexErrMsgTxt("Invalid output matrix size.");
+    }
+    ndims = mxGetM(N_IN)*mxGetN(N_IN);
+    if (ndims<2 || ndims>3) {
+      mexErrMsgTxt("Output matrix can only be 2-D or 3-D.");
+    }
+    if (!mxIsDouble(SZ_IN) || mxIsComplex(SZ_IN) || mxGetNumberOfDimensions(SZ_IN)>2 || mxGetM(SZ_IN)*mxGetN(SZ_IN)!=ndims) {
+      mexErrMsgTxt("Invalid block size.");
+    }
+    if (nrhs == 4) {
+      if (!mxIsDouble(S_IN) || mxIsComplex(S_IN) || mxGetNumberOfDimensions(S_IN)>2 || mxGetM(S_IN)*mxGetN(S_IN)!=ndims) {
+        mexErrMsgTxt("Invalid step size.");
+      }
+    }
+
+
+    /* Get parameters */
+
+    s = mxGetPr(N_IN);
+    if (s[0]<1 || s[1]<1 || (ndims==3 && s[2]<1)) {
+      mexErrMsgTxt("Invalid output matrix size.");
+    }
+    n[0] = (mwSize)(s[0] + 0.01);
+    n[1] = (mwSize)(s[1] + 0.01);
+    n[2] = ndims==3 ? (mwSize)(s[2] + 0.01) : 1;
+
+    s = mxGetPr(SZ_IN);
+    if (s[0]<1 || s[1]<1 || (ndims==3 && s[2]<1)) {
+      mexErrMsgTxt("Invalid block size.");
+    }
+    sz[0] = (mwSize)(s[0] + 0.01);
+    sz[1] = (mwSize)(s[1] + 0.01);
+    sz[2] = ndims==3 ? (mwSize)(s[2] + 0.01) : 1;
+
+    if (nrhs == 4) {
+      s = mxGetPr(S_IN);
+      if (s[0]<1 || s[1]<1 || (ndims==3 && s[2]<1)) {
+        mexErrMsgTxt("Invalid step size.");
+      }
+      stepsize[0] = (mwSize)(s[0] + 0.01);
+      stepsize[1] = (mwSize)(s[1] + 0.01);
+      stepsize[2] = ndims==3 ? (mwSize)(s[2] + 0.01) : 1;
+    }
+    else {
+      stepsize[0] = stepsize[1] = stepsize[2] = 1;
+    }
+
+    if (n[0]<sz[0] || n[1]<sz[1] || (ndims==3 && n[2]<sz[2])) {
+      mexErrMsgTxt("Block size too large.");
+    }
+
+    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)) {
+      mexErrMsgTxt("Invalid number of columns in B. Please use IM2COLSTEP to compute B.");
+    }
+
+
+    /* Create a matrix for the return argument */
+
+    X_OUT = mxCreateNumericArray(ndims, n, mxDOUBLE_CLASS, mxREAL);
+
+
+    /* Assign pointers */
+
+    b = mxGetPr(B_IN);
+    x = mxGetPr(X_OUT);
+
+
+    /* Do the actual computation */
+
+    blocknum = 0;
+
+    /* iterate over all blocks */
+    for (k=0; k<=n[2]-sz[2]; k+=stepsize[2]) {
+      for (j=0; j<=n[1]-sz[1]; j+=stepsize[1]) {
+        for (i=0; i<=n[0]-sz[0]; i+=stepsize[0]) {
+
+          /* add single block */
+          for (m=0; m<sz[2]; m++) {
+            for (l=0; l<sz[1]; l++) {
+              for (t=0; t<sz[0]; t++) {
+                (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];
+              }
+            }
+          }
+          blocknum++;
+
+        }
+      }
+    }
+
+    return;
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Problems/private/col2imstep.m	Thu Mar 25 14:05:27 2010 +0000
@@ -0,0 +1,25 @@
+%COL2IMSTEP Rearrange matrix columns into blocks.
+%  A = COL2IMSTEP(B,[MM NN],[N1 N2]) rearranges the columns of B into
+%  sliding N1-by-N2 blocks producing the matrix A of size MM-by-NN. B is
+%  usually the result of calling IM2COLSTEP(...) with a stepsize of 1, or
+%  using Matlab's IM2COL(..,'sliding'). Overlapping blocks are summed in A.
+%
+%  A = COL2IMSTEP(B,[MM NN],[N1 N2],[S1 S2]) arranges the blocks in A with
+%  a step size of (S1,S2) between them. The first block is at A(1:N1,1:N2),
+%  and the rest are at A((1:N1)+i*S1,(1:N2)+j*S2). Overlapping blocks are
+%  summed in A. Note that B is usually the result of calling
+%  IM2COLSTEP(...) with a stepsize of [S1 S2].
+%
+%  A = IM2COLSTEP(B,[MM NN KK],[N1 N2 N3],[S1 S2 S3]) generates a 3-D
+%  output matrix A. The step size [S1 S2 S3] may be ommitted, and defaults
+%  to [1 1 1].
+%
+%  See also IM2COLSTEP, IM2COL, COUNTCOVER.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  August 2009
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Problems/private/im2colstep.c	Thu Mar 25 14:05:27 2010 +0000
@@ -0,0 +1,131 @@
+/**************************************************************************
+ *
+ * File name: im2colstep.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 31.8.2009
+ *
+ *************************************************************************/
+
+
+#include "mex.h"
+#include <string.h>
+
+
+/* Input Arguments */
+
+#define	X_IN	 prhs[0]
+#define SZ_IN  prhs[1]
+#define S_IN   prhs[2]
+
+
+/* Output Arguments */
+
+#define	B_OUT	plhs[0]
+
+
+void mexFunction(int nlhs, mxArray *plhs[],
+		             int nrhs, const mxArray*prhs[])
+
+{
+    double *x, *b, *s;
+    mwSize sz[3], stepsize[3], n[3], ndims;
+    mwIndex i, j, k, l, m, blocknum;
+
+
+    /* Check for proper number of arguments */
+
+    if (nrhs < 2 || nrhs > 3) {
+      mexErrMsgTxt("Invalid number of input arguments.");
+    } else if (nlhs > 1) {
+      mexErrMsgTxt("Too many output arguments.");
+    }
+
+
+    /* Check the the input dimensions */
+
+    ndims = mxGetNumberOfDimensions(X_IN);
+
+    if (!mxIsDouble(X_IN) || mxIsComplex(X_IN) || ndims>3) {
+      mexErrMsgTxt("X should be a 2-D or 3-D double matrix.");
+    }
+    if (!mxIsDouble(SZ_IN) || mxIsComplex(SZ_IN) || mxGetNumberOfDimensions(SZ_IN)>2 || mxGetM(SZ_IN)*mxGetN(SZ_IN)!=ndims) {
+      mexErrMsgTxt("Invalid block size.");
+    }
+    if (nrhs == 3) {
+      if (!mxIsDouble(S_IN) || mxIsComplex(S_IN) || mxGetNumberOfDimensions(S_IN)>2 || mxGetM(S_IN)*mxGetN(S_IN)!=ndims) {
+        mexErrMsgTxt("Invalid step size.");
+      }
+    }
+
+
+    /* Get parameters */
+
+    s = mxGetPr(SZ_IN);
+    if (s[0]<1 || s[1]<1 || (ndims==3 && s[2]<1)) {
+      mexErrMsgTxt("Invalid block size.");
+    }
+    sz[0] = (mwSize)(s[0] + 0.01);
+    sz[1] = (mwSize)(s[1] + 0.01);
+    sz[2] = ndims==3 ? (mwSize)(s[2] + 0.01) : 1;
+
+    if (nrhs == 3) {
+      s = mxGetPr(S_IN);
+      if (s[0]<1 || s[1]<1 || (ndims==3 && s[2]<1)) {
+        mexErrMsgTxt("Invalid step size.");
+      }
+      stepsize[0] = (mwSize)(s[0] + 0.01);
+      stepsize[1] = (mwSize)(s[1] + 0.01);
+      stepsize[2] = ndims==3 ? (mwSize)(s[2] + 0.01) : 1;
+    }
+    else {
+      stepsize[0] = stepsize[1] = stepsize[2] = 1;
+    }
+
+    n[0] = (mxGetDimensions(X_IN))[0];
+    n[1] = (mxGetDimensions(X_IN))[1];
+    n[2] = ndims==3 ? (mxGetDimensions(X_IN))[2] : 1;
+
+    if (n[0]<sz[0] || n[1]<sz[1] || (ndims==3 && n[2]<sz[2])) {
+      mexErrMsgTxt("Block size too large.");
+    }
+
+
+    /* Create a matrix for the return argument */
+
+    B_OUT = mxCreateDoubleMatrix(sz[0]*sz[1]*sz[2], ((n[0]-sz[0])/stepsize[0]+1)*((n[1]-sz[1])/stepsize[1]+1)*((n[2]-sz[2])/stepsize[2]+1), mxREAL);
+
+
+    /* Assign pointers */
+
+    x = mxGetPr(X_IN);
+    b = mxGetPr(B_OUT);
+
+
+    /* Do the actual computation */
+
+    blocknum = 0;
+
+    /* iterate over all blocks */
+    for (k=0; k<=n[2]-sz[2]; k+=stepsize[2]) {
+      for (j=0; j<=n[1]-sz[1]; j+=stepsize[1]) {
+        for (i=0; i<=n[0]-sz[0]; i+=stepsize[0]) {
+
+          /* copy single block */
+          for (m=0; m<sz[2]; m++) {
+            for (l=0; l<sz[1]; l++) {
+              memcpy(b + blocknum*sz[0]*sz[1]*sz[2] + m*sz[0]*sz[1] + l*sz[0], x+(k+m)*n[0]*n[1]+(j+l)*n[0]+i, sz[0]*sizeof(double));
+            }
+          }
+          blocknum++;
+
+        }
+      }
+    }
+
+    return;
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Problems/private/im2colstep.m	Thu Mar 25 14:05:27 2010 +0000
@@ -0,0 +1,31 @@
+%IM2COLSTEP Rearrange matrix blocks into columns.
+%  B = IM2COLSTEP(A,[N1 N2]) converts each sliding N1-by-N2 block of the
+%  2-D matrix A into a column of B, with no zero padding. B has N1*N2 rows
+%  and will contain as many columns as there are N1-by-N2 neighborhoods in
+%  A. Each column of B contains a neighborhood of A reshaped as NHOOD(:),
+%  where NHOOD is a matrix containing an N1-by-N2 neighborhood of A.
+%
+%  B = IM2COLSTEP(A,[N1 N2],[S1 S2]) extracts neighborhoods of A with a
+%  step size of (S1,S2) between them. The first extracted neighborhood is
+%  A(1:N1,1:N2), and the rest are of the form A((1:N1)+i*S1,(1:N2)+j*S2).
+%  Note that to ensure coverage of all A by neighborhoods,
+%  (size(A,i)-Ni)/Si must be whole for i=1,2. The default function behavior
+%  corresponds to [S1 S2] = [1 1]. Setting S1>=N1 and S2>=N2 results in no
+%  overlap between the neighborhoods.
+%
+%  B = IM2COLSTEP(A,[N1 N2 N3],[S1 S2 S3]) operates on a 3-D matrix A. The
+%  step size [S1 S2 S3] may be ommitted, and defaults to [1 1 1].
+%
+%  Note: the call IM2COLSTEP(A,[N1 N2]) produces the same output as
+%  Matlab's IM2COL(A,[N1 N2],'sliding'). However, it is significantly
+%  faster.
+%
+%  See also COL2IMSTEP, IM2COL, COUNTCOVER.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  August 2009
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Problems/private/mexutils.c	Thu Mar 25 14:05:27 2010 +0000
@@ -0,0 +1,79 @@
+/**************************************************************************
+ *
+ * File name: mexutils.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 15.8.2009
+ *
+ *************************************************************************/
+
+#include "mexutils.h"
+#include <math.h>
+
+
+
+/* verify that the mxArray contains a double matrix */
+
+void checkmatrix(const mxArray *param, char *fname, char *pname)
+{
+  char errmsg[100];
+  sprintf(errmsg, "%.15s requires that %.25s be a double matrix.", fname, pname);
+  if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2) {
+    mexErrMsgTxt(errmsg);
+  }
+}
+
+
+/* verify that the mxArray contains a 1-D double vector */
+
+void checkvector(const mxArray *param, char *fname, char *pname)
+{
+  char errmsg[100];
+  sprintf(errmsg, "%.15s requires that %.25s be a double vector.", fname, pname);
+  if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2 || (mxGetM(param)!=1 && mxGetN(param)!=1)) {
+    mexErrMsgTxt(errmsg);
+  }
+}
+
+
+/* verify that the mxArray contains a double scalar */
+
+void checkscalar(const mxArray *param, char *fname, char *pname)
+{
+  char errmsg[100];
+  sprintf(errmsg, "%.15s requires that %.25s be a double scalar.", fname, pname);
+  if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2 ||
+      mxGetM(param)!=1 || mxGetN(param)!=1)
+  {
+    mexErrMsgTxt(errmsg);
+  }
+}
+
+
+/* verify that the mxArray contains a sparse matrix */
+
+void checksparse(const mxArray *param, char *fname, char *pname)
+{
+  char errmsg[100];
+  sprintf(errmsg, "%.15s requires that %.25s be sparse.", fname, pname);
+  if (!mxIsSparse(param)) {
+    mexErrMsgTxt(errmsg);
+  }
+}
+
+
+/* verify that the mxArray contains a 1-dimensional cell array */
+
+void checkcell_1d(const mxArray *param, char *fname, char *pname)
+{
+  char errmsg[100];
+  sprintf(errmsg, "%.15s requires that %.25s be a 1-D cell array.", fname, pname);
+  if (!mxIsCell(param) || mxGetNumberOfDimensions(param)>2 || (mxGetM(param)!=1 && mxGetN(param)!=1)) {
+    mexErrMsgTxt(errmsg);
+  }
+}
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Problems/private/mexutils.h	Thu Mar 25 14:05:27 2010 +0000
@@ -0,0 +1,103 @@
+/**************************************************************************
+ *
+ * File name: mexutils.h
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 18.8.2009
+ *
+ * Utility functions for MEX files.
+ *
+ *************************************************************************/
+
+
+#ifndef __MEX_UTILS_H__
+#define __MEX_UTILS_H__
+
+#include "mex.h"
+
+
+
+/**************************************************************************
+ * Function checkmatrix:
+ *
+ * Verify that the specified mxArray is real, of type double, and has
+ * no more than two dimensions. If not, an error message is printed
+ * and the mex file terminates.
+ *
+ * Parameters:
+ *   param - the mxArray to be checked
+ *   fname - the name of the function where the error occured (15 characters or less)
+ *   pname - the name of the parameter (25 characters or less)
+ *
+ **************************************************************************/
+void checkmatrix(const mxArray *param, char *fname, char *pname);
+
+
+/**************************************************************************
+ * Function checkvector:
+ *
+ * Verify that the specified mxArray is 1-D, real, and of type double. The
+ * vector may be a column or row vector. Otherwise, an error message is
+ * printed and the mex file terminates.
+ *
+ * Parameters:
+ *   param - the mxArray to be checked
+ *   fname - the name of the function where the error occured (15 characters or less)
+ *   pname - the name of the parameter (25 characters or less)
+ *
+ **************************************************************************/
+void checkvector(const mxArray *param, char *fname, char *pname);
+
+
+/**************************************************************************
+ * Function checkscalar:
+ *
+ * Verify that the specified mxArray represents a real double scalar value.
+ * If not, an error message is printed and the mex file terminates.
+ *
+ * Parameters:
+ *   param - the mxArray to be checked
+ *   fname - the name of the function where the error occured (15 characters or less)
+ *   pname - the name of the parameter (25 characters or less)
+ *
+ **************************************************************************/
+void checkscalar(const mxArray *param, char *fname, char *pname);
+
+
+/**************************************************************************
+ * Function checksparse:
+ *
+ * Verify that the specified mxArray contains a sparse matrix. If not,
+ * an error message is printed and the mex file terminates.
+ *
+ * Parameters:
+ *   param - the mxArray to be checked
+ *   fname - the name of the function where the error occured (15 characters or less)
+ *   pname - the name of the parameter (25 characters or less)
+ *
+ **************************************************************************/
+void checksparse(const mxArray *param, char *fname, char *pname);
+
+
+/**************************************************************************
+ * Function checkcell_1d:
+ *
+ * Verify that the specified mxArray is a 1-D cell array. The cell array
+ * may be arranged as either a column or a row. If not, an error message
+ * is printed and the mex file terminates.
+ *
+ * Parameters:
+ *   param - the mxArray to be checked
+ *   fname - the name of the function where the error occured (15 characters or less)
+ *   pname - the name of the parameter (25 characters or less)
+ *
+ **************************************************************************/
+void checkcell_1d(const mxArray *param, char *fname, char *pname);
+
+
+#endif
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Problems/private/sprow.c	Thu Mar 25 14:05:27 2010 +0000
@@ -0,0 +1,94 @@
+/**************************************************************************
+ *
+ * File name: sprow.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 24.8.2009
+ *
+ *************************************************************************/
+
+
+#include "mex.h"
+#include "mexutils.h"
+
+
+/* Input Arguments */
+
+#define	A_IN	prhs[0]
+#define J_IN  prhs[1]
+
+
+/* Output Arguments */
+
+#define	X_OUT	  plhs[0]
+#define	ID_OUT	plhs[1]
+
+
+void mexFunction(int nlhs, mxArray *plhs[],
+		             int nrhs, const mxArray*prhs[])
+
+{
+  double *pr, *x, *id, rowid;
+  mwIndex *ir, *jc;
+  mwSize m, n;
+  mwIndex i, j, k, l, rowlen;
+
+  if (nrhs != 2) {
+    mexErrMsgTxt("GETSPROW requires two input arguments.");
+  } else if (nlhs > 2) {
+    mexErrMsgTxt("Too many output arguments.");
+  }
+
+  checkmatrix(A_IN, "GETSPROW", "A");
+  checksparse(A_IN, "GETSPROW", "A");
+  checkscalar(J_IN, "GETSPROW", "J");
+
+  m = mxGetM(A_IN);
+  n = mxGetN(A_IN);
+
+  rowid = mxGetScalar(J_IN);
+  if (rowid < 0) {
+    mexErrMsgTxt("Invalid row index.");
+  }
+  j = (mwIndex)(rowid + 1e-2);
+  if (j<1 || j>m) {
+    mexErrMsgTxt("Row index is out of range.");
+  }
+  j--;
+
+  pr = mxGetPr(A_IN);
+  ir = mxGetIr(A_IN);
+  jc = mxGetJc(A_IN);
+
+  /* Determine length of row */
+  rowlen = 0;
+  for (i=0; i<jc[n]; ++i) {
+    rowlen += (ir[i]==j);
+  }
+
+  /* Allocate output parameters */
+  X_OUT = mxCreateDoubleMatrix(1, rowlen, mxREAL);
+  ID_OUT = mxCreateDoubleMatrix(1, rowlen, mxREAL);
+
+  x = mxGetPr(X_OUT);
+  id = mxGetPr(ID_OUT);
+
+  /* Compute j-th row */
+  k=0;
+  for (l=1; l<=n; ++l) {
+    i = jc[l-1];
+    while (i<jc[l] && ir[i]<j) {
+      i++;
+    }
+    if (i<jc[l] && ir[i]==j) {
+      x[k] = pr[i];
+      id[k] = l;
+      k++;
+    }
+  }
+
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Problems/private/sprow.m	Thu Mar 25 14:05:27 2010 +0000
@@ -0,0 +1,19 @@
+%SPROW Extract row of sparse matrix.
+%  X = SPROW(A,J) where A is a sparse matrix, returns the nonzero values in
+%  the row A(J,:).
+%
+%  [X,ID] = SPROW(A,J) also returns the column indices of the nonzeros.
+%
+%  Note that the call [X,ID] = SPROW(A,J) is equivalent (but more efficient
+%  than) the Matlab code
+%
+%    IDS = find(A(J,:));
+%    X = A(J,IDS);
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  August 2009