1 | #include "mex.h" |
---|
2 | |
---|
3 | /* Returns the interpolated output value at index, i.e. (index, value) |
---|
4 | * using the known (inp, outp) pairs. |
---|
5 | */ |
---|
6 | void table1(double *inp, double *outp, int length, double index, |
---|
7 | int output_length, double *value_return, int stride) |
---|
8 | { |
---|
9 | int i,j; |
---|
10 | double k; |
---|
11 | /* find the first element of inp greater than index */ |
---|
12 | #if 1 |
---|
13 | for(i=0;i<length;i++) { |
---|
14 | if(inp[i] > index) break; |
---|
15 | } |
---|
16 | #else |
---|
17 | /* binary search */ |
---|
18 | int low = 0, high = length; |
---|
19 | while(low < high) { |
---|
20 | i = (high+low)/2; |
---|
21 | if(inp[i] >= index) high = i; |
---|
22 | else if(inp[i] < index) low = i+1; |
---|
23 | } |
---|
24 | #endif |
---|
25 | /* too small or large? guess the extreme value */ |
---|
26 | if(i == 0) { |
---|
27 | k = 1; |
---|
28 | for(j=0;j<output_length;j++) { |
---|
29 | *value_return = outp[i]; |
---|
30 | value_return += stride; |
---|
31 | outp += length; |
---|
32 | } |
---|
33 | } |
---|
34 | else if(i == length) { |
---|
35 | k = 0; |
---|
36 | for(j=0;j<output_length;j++) { |
---|
37 | *value_return = outp[i-1]; |
---|
38 | value_return += stride; |
---|
39 | outp += length; |
---|
40 | } |
---|
41 | } |
---|
42 | else { |
---|
43 | /* interpolate between this element and the previous one */ |
---|
44 | k = (index - inp[i-1]) / (inp[i] - inp[i-1]); |
---|
45 | for(j=0;j<output_length;j++) { |
---|
46 | *value_return = k*outp[i] + (1-k)*outp[i-1]; |
---|
47 | value_return += stride; |
---|
48 | outp += length; |
---|
49 | } |
---|
50 | } |
---|
51 | } |
---|
52 | |
---|
53 | void mexFunction(int nlhs, mxArray *plhs[], |
---|
54 | int nrhs, const mxArray *prhs[]) |
---|
55 | { |
---|
56 | double *pr; |
---|
57 | int length, data_length, nx, i; |
---|
58 | |
---|
59 | if((nrhs != 2) || (nlhs > 1)) mexErrMsgTxt("Usage: y = table1(tab, x)"); |
---|
60 | pr = mxGetPr(prhs[0]); |
---|
61 | length = mxGetM(prhs[0]); |
---|
62 | data_length = mxGetN(prhs[0])-1; |
---|
63 | nx = mxGetM(prhs[1]); |
---|
64 | if(mxGetN(prhs[1]) > nx) nx = mxGetN(prhs[1]); |
---|
65 | plhs[0] = mxCreateDoubleMatrix(nx, data_length, REAL); |
---|
66 | for(i=0;i<nx;i++) { |
---|
67 | table1(pr, pr+length, length, mxGetPr(prhs[1])[i], |
---|
68 | data_length, mxGetPr(plhs[0])+i, nx); |
---|
69 | } |
---|
70 | } |
---|