matthiasm@8: /* This is based on
matthiasm@8: http://www.mathworks.com/access/helpdesk/help/techdoc/matlab_external/ch04cr12.shtml
matthiasm@8: 
matthiasm@8: See rectintSparse.m for the matlab version of this code.
matthiasm@8: 
matthiasm@8: */
matthiasm@8: 
matthiasm@8: #include <math.h> /* Needed for the ceil() prototype. */
matthiasm@8: #include "mex.h"
matthiasm@8: #include <stdio.h>
matthiasm@8: 
matthiasm@8: /* If you are using a compiler that equates NaN to be zero, you 
matthiasm@8:  * must compile this example using the flag  -DNAN_EQUALS_ZERO.
matthiasm@8:  * For example:
matthiasm@8:  *
matthiasm@8:  *     mex -DNAN_EQUALS_ZERO fulltosparse.c
matthiasm@8:  *
matthiasm@8:  * This will correctly define the IsNonZero macro for your C
matthiasm@8:  * compiler.
matthiasm@8:  */
matthiasm@8: 
matthiasm@8: #if defined(NAN_EQUALS_ZERO)
matthiasm@8: #define IsNonZero(d) ((d) != 0.0 || mxIsNaN(d))
matthiasm@8: #else
matthiasm@8: #define IsNonZero(d) ((d) != 0.0)
matthiasm@8: #endif
matthiasm@8: 
matthiasm@8: #define MAX(x,y)  ((x)>(y) ? (x) : (y))
matthiasm@8: #define MIN(x,y)  ((x)<(y) ? (x) : (y))
matthiasm@8: 
matthiasm@8: void mexFunction(
matthiasm@8:         int nlhs,       mxArray *plhs[],
matthiasm@8:         int nrhs, const mxArray *prhs[]
matthiasm@8:         )
matthiasm@8: {
matthiasm@8:   /* Declare variables. */
matthiasm@8:   int j,k,m,n,nzmax,*irs,*jcs, *irs2, *jcs2;
matthiasm@8:   double *overlap, *overlap2, tmp, areaA, areaB;
matthiasm@8:   double percent_sparse;
matthiasm@8:   double *leftA, *rightA, *topA, *bottomA;
matthiasm@8:   double *leftB, *rightB, *topB, *bottomB;
matthiasm@8:   double *verbose;
matthiasm@8: 
matthiasm@8:   /* Get the size and pointers to input data. */
matthiasm@8:   m = MAX(mxGetM(prhs[0]), mxGetN(prhs[0]));
matthiasm@8:   n = MAX(mxGetM(prhs[4]), mxGetN(prhs[4]));
matthiasm@8:   /* printf("A=%d, B=%d\n", m, n); */
matthiasm@8: 
matthiasm@8:   leftA = mxGetPr(prhs[0]);
matthiasm@8:   rightA = mxGetPr(prhs[1]);
matthiasm@8:   topA = mxGetPr(prhs[2]);
matthiasm@8:   bottomA = mxGetPr(prhs[3]);
matthiasm@8: 
matthiasm@8:   leftB = mxGetPr(prhs[4]);
matthiasm@8:   rightB = mxGetPr(prhs[5]);
matthiasm@8:   topB = mxGetPr(prhs[6]);
matthiasm@8:   bottomB = mxGetPr(prhs[7]);
matthiasm@8: 
matthiasm@8:   verbose = mxGetPr(prhs[8]);
matthiasm@8: 
matthiasm@8:     /* Allocate space for sparse matrix. 
matthiasm@8:      * NOTE:  Assume at most 20% of the data is sparse.  Use ceil
matthiasm@8:      * to cause it to round up. 
matthiasm@8:      */
matthiasm@8: 
matthiasm@8:   percent_sparse = 0.01;
matthiasm@8:   nzmax = (int)ceil((double)m*(double)n*percent_sparse);
matthiasm@8: 
matthiasm@8:   plhs[0] = mxCreateSparse(m,n,nzmax,0);
matthiasm@8:   overlap  = mxGetPr(plhs[0]);
matthiasm@8:   irs = mxGetIr(plhs[0]);
matthiasm@8:   jcs = mxGetJc(plhs[0]);
matthiasm@8: 
matthiasm@8:   plhs[1] = mxCreateSparse(m,n,nzmax,0);
matthiasm@8:   overlap2  = mxGetPr(plhs[1]);
matthiasm@8:   irs2 = mxGetIr(plhs[1]);
matthiasm@8:   jcs2 = mxGetJc(plhs[1]);
matthiasm@8: 
matthiasm@8:     
matthiasm@8:   /* Assign nonzeros. */
matthiasm@8:   k = 0; 
matthiasm@8:   for (j = 0; (j < n); j++) {
matthiasm@8:     int i;
matthiasm@8:     jcs[j] = k; 
matthiasm@8:     jcs2[j] = k; 
matthiasm@8:     for (i = 0; (i < m); i++) {
matthiasm@8:       tmp = (MAX(0, MIN(rightA[i], rightB[j]) - MAX(leftA[i], leftB[j]) )) * 
matthiasm@8: 	(MAX(0, MIN(topA[i], topB[j]) - MAX(bottomA[i], bottomB[j]) ));
matthiasm@8:       
matthiasm@8:       if (*verbose) {
matthiasm@8: 	printf("j=%d,i=%d,tmp=%5.3f\n", j,i,tmp);
matthiasm@8:       }
matthiasm@8: 
matthiasm@8:       if (IsNonZero(tmp)) {
matthiasm@8: 
matthiasm@8:         /* Check to see if non-zero element will fit in 
matthiasm@8:          * allocated output array.  If not, increase
matthiasm@8:          * percent_sparse by 20%, recalculate nzmax, and augment
matthiasm@8:          * the sparse array.
matthiasm@8:          */
matthiasm@8:         if (k >= nzmax) {
matthiasm@8:           int oldnzmax = nzmax;
matthiasm@8:           percent_sparse += 0.2;
matthiasm@8:           nzmax = (int)ceil((double)m*(double)n*percent_sparse);
matthiasm@8: 
matthiasm@8:           /* Make sure nzmax increases atleast by 1. */
matthiasm@8:           if (oldnzmax == nzmax) 
matthiasm@8:             nzmax++;
matthiasm@8: 	  printf("reallocating from %d to %d\n", oldnzmax, nzmax);
matthiasm@8: 
matthiasm@8:           mxSetNzmax(plhs[0], nzmax); 
matthiasm@8:           mxSetPr(plhs[0], mxRealloc(overlap, nzmax*sizeof(double)));
matthiasm@8:           mxSetIr(plhs[0], mxRealloc(irs, nzmax*sizeof(int)));
matthiasm@8:           overlap  = mxGetPr(plhs[0]);
matthiasm@8:           irs = mxGetIr(plhs[0]);
matthiasm@8: 
matthiasm@8:           mxSetNzmax(plhs[1], nzmax); 
matthiasm@8:           mxSetPr(plhs[1], mxRealloc(overlap2, nzmax*sizeof(double)));
matthiasm@8:           mxSetIr(plhs[1], mxRealloc(irs2, nzmax*sizeof(int)));
matthiasm@8:           overlap2  = mxGetPr(plhs[1]);
matthiasm@8:           irs2 = mxGetIr(plhs[1]);
matthiasm@8:         }
matthiasm@8: 
matthiasm@8:         overlap[k] = tmp;
matthiasm@8:         irs[k] = i;
matthiasm@8: 	
matthiasm@8: 	areaA = (rightA[i]-leftA[i])*(topA[i]-bottomA[i]);
matthiasm@8: 	areaB = (rightB[j]-leftB[j])*(topB[j]-bottomB[j]);
matthiasm@8: 	overlap2[k] = MIN(tmp/areaA, tmp/areaB);
matthiasm@8: 	irs2[k] = i;
matthiasm@8: 
matthiasm@8:         k++;
matthiasm@8:       } /* IsNonZero */
matthiasm@8:     } /* for i */
matthiasm@8:   }
matthiasm@8:   jcs[n] = k;
matthiasm@8:   jcs2[n] = k;
matthiasm@8:   
matthiasm@8: }
matthiasm@8: 
matthiasm@8: 
matthiasm@8: 
matthiasm@8: 
matthiasm@8: 
matthiasm@8: 
matthiasm@8: 
matthiasm@8: