Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/KPMtools/rectintSparseLoopC.c @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 /* This is based on
wolffd@0 2 http://www.mathworks.com/access/helpdesk/help/techdoc/matlab_external/ch04cr12.shtml
wolffd@0 3
wolffd@0 4 See rectintSparse.m for the matlab version of this code.
wolffd@0 5
wolffd@0 6 */
wolffd@0 7
wolffd@0 8 #include <math.h> /* Needed for the ceil() prototype. */
wolffd@0 9 #include "mex.h"
wolffd@0 10 #include <stdio.h>
wolffd@0 11
wolffd@0 12 /* If you are using a compiler that equates NaN to be zero, you
wolffd@0 13 * must compile this example using the flag -DNAN_EQUALS_ZERO.
wolffd@0 14 * For example:
wolffd@0 15 *
wolffd@0 16 * mex -DNAN_EQUALS_ZERO fulltosparse.c
wolffd@0 17 *
wolffd@0 18 * This will correctly define the IsNonZero macro for your C
wolffd@0 19 * compiler.
wolffd@0 20 */
wolffd@0 21
wolffd@0 22 #if defined(NAN_EQUALS_ZERO)
wolffd@0 23 #define IsNonZero(d) ((d) != 0.0 || mxIsNaN(d))
wolffd@0 24 #else
wolffd@0 25 #define IsNonZero(d) ((d) != 0.0)
wolffd@0 26 #endif
wolffd@0 27
wolffd@0 28 #define MAX(x,y) ((x)>(y) ? (x) : (y))
wolffd@0 29 #define MIN(x,y) ((x)<(y) ? (x) : (y))
wolffd@0 30
wolffd@0 31 void mexFunction(
wolffd@0 32 int nlhs, mxArray *plhs[],
wolffd@0 33 int nrhs, const mxArray *prhs[]
wolffd@0 34 )
wolffd@0 35 {
wolffd@0 36 /* Declare variables. */
wolffd@0 37 int j,k,m,n,nzmax,*irs,*jcs, *irs2, *jcs2;
wolffd@0 38 double *overlap, *overlap2, tmp, areaA, areaB;
wolffd@0 39 double percent_sparse;
wolffd@0 40 double *leftA, *rightA, *topA, *bottomA;
wolffd@0 41 double *leftB, *rightB, *topB, *bottomB;
wolffd@0 42 double *verbose;
wolffd@0 43
wolffd@0 44 /* Get the size and pointers to input data. */
wolffd@0 45 m = MAX(mxGetM(prhs[0]), mxGetN(prhs[0]));
wolffd@0 46 n = MAX(mxGetM(prhs[4]), mxGetN(prhs[4]));
wolffd@0 47 /* printf("A=%d, B=%d\n", m, n); */
wolffd@0 48
wolffd@0 49 leftA = mxGetPr(prhs[0]);
wolffd@0 50 rightA = mxGetPr(prhs[1]);
wolffd@0 51 topA = mxGetPr(prhs[2]);
wolffd@0 52 bottomA = mxGetPr(prhs[3]);
wolffd@0 53
wolffd@0 54 leftB = mxGetPr(prhs[4]);
wolffd@0 55 rightB = mxGetPr(prhs[5]);
wolffd@0 56 topB = mxGetPr(prhs[6]);
wolffd@0 57 bottomB = mxGetPr(prhs[7]);
wolffd@0 58
wolffd@0 59 verbose = mxGetPr(prhs[8]);
wolffd@0 60
wolffd@0 61 /* Allocate space for sparse matrix.
wolffd@0 62 * NOTE: Assume at most 20% of the data is sparse. Use ceil
wolffd@0 63 * to cause it to round up.
wolffd@0 64 */
wolffd@0 65
wolffd@0 66 percent_sparse = 0.01;
wolffd@0 67 nzmax = (int)ceil((double)m*(double)n*percent_sparse);
wolffd@0 68
wolffd@0 69 plhs[0] = mxCreateSparse(m,n,nzmax,0);
wolffd@0 70 overlap = mxGetPr(plhs[0]);
wolffd@0 71 irs = mxGetIr(plhs[0]);
wolffd@0 72 jcs = mxGetJc(plhs[0]);
wolffd@0 73
wolffd@0 74 plhs[1] = mxCreateSparse(m,n,nzmax,0);
wolffd@0 75 overlap2 = mxGetPr(plhs[1]);
wolffd@0 76 irs2 = mxGetIr(plhs[1]);
wolffd@0 77 jcs2 = mxGetJc(plhs[1]);
wolffd@0 78
wolffd@0 79
wolffd@0 80 /* Assign nonzeros. */
wolffd@0 81 k = 0;
wolffd@0 82 for (j = 0; (j < n); j++) {
wolffd@0 83 int i;
wolffd@0 84 jcs[j] = k;
wolffd@0 85 jcs2[j] = k;
wolffd@0 86 for (i = 0; (i < m); i++) {
wolffd@0 87 tmp = (MAX(0, MIN(rightA[i], rightB[j]) - MAX(leftA[i], leftB[j]) )) *
wolffd@0 88 (MAX(0, MIN(topA[i], topB[j]) - MAX(bottomA[i], bottomB[j]) ));
wolffd@0 89
wolffd@0 90 if (*verbose) {
wolffd@0 91 printf("j=%d,i=%d,tmp=%5.3f\n", j,i,tmp);
wolffd@0 92 }
wolffd@0 93
wolffd@0 94 if (IsNonZero(tmp)) {
wolffd@0 95
wolffd@0 96 /* Check to see if non-zero element will fit in
wolffd@0 97 * allocated output array. If not, increase
wolffd@0 98 * percent_sparse by 20%, recalculate nzmax, and augment
wolffd@0 99 * the sparse array.
wolffd@0 100 */
wolffd@0 101 if (k >= nzmax) {
wolffd@0 102 int oldnzmax = nzmax;
wolffd@0 103 percent_sparse += 0.2;
wolffd@0 104 nzmax = (int)ceil((double)m*(double)n*percent_sparse);
wolffd@0 105
wolffd@0 106 /* Make sure nzmax increases atleast by 1. */
wolffd@0 107 if (oldnzmax == nzmax)
wolffd@0 108 nzmax++;
wolffd@0 109 printf("reallocating from %d to %d\n", oldnzmax, nzmax);
wolffd@0 110
wolffd@0 111 mxSetNzmax(plhs[0], nzmax);
wolffd@0 112 mxSetPr(plhs[0], mxRealloc(overlap, nzmax*sizeof(double)));
wolffd@0 113 mxSetIr(plhs[0], mxRealloc(irs, nzmax*sizeof(int)));
wolffd@0 114 overlap = mxGetPr(plhs[0]);
wolffd@0 115 irs = mxGetIr(plhs[0]);
wolffd@0 116
wolffd@0 117 mxSetNzmax(plhs[1], nzmax);
wolffd@0 118 mxSetPr(plhs[1], mxRealloc(overlap2, nzmax*sizeof(double)));
wolffd@0 119 mxSetIr(plhs[1], mxRealloc(irs2, nzmax*sizeof(int)));
wolffd@0 120 overlap2 = mxGetPr(plhs[1]);
wolffd@0 121 irs2 = mxGetIr(plhs[1]);
wolffd@0 122 }
wolffd@0 123
wolffd@0 124 overlap[k] = tmp;
wolffd@0 125 irs[k] = i;
wolffd@0 126
wolffd@0 127 areaA = (rightA[i]-leftA[i])*(topA[i]-bottomA[i]);
wolffd@0 128 areaB = (rightB[j]-leftB[j])*(topB[j]-bottomB[j]);
wolffd@0 129 overlap2[k] = MIN(tmp/areaA, tmp/areaB);
wolffd@0 130 irs2[k] = i;
wolffd@0 131
wolffd@0 132 k++;
wolffd@0 133 } /* IsNonZero */
wolffd@0 134 } /* for i */
wolffd@0 135 }
wolffd@0 136 jcs[n] = k;
wolffd@0 137 jcs2[n] = k;
wolffd@0 138
wolffd@0 139 }
wolffd@0 140
wolffd@0 141
wolffd@0 142
wolffd@0 143
wolffd@0 144
wolffd@0 145
wolffd@0 146
wolffd@0 147