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
|