Implementing a tensor product formula

I would like to use C to implement the following formula:

$$\mathbf S(u,v)=\sum_{r=i-p}^{i}\sum_{s=j-q}^{j}N_{r,p}(u)N_{s,q}(v)\mathbf P_{r,s}$$

Namely,

$$
\left(N_{i-p,p}(u),\cdots,N_{i,p}(u)\right)_{1\times(p+1)}
\begin{pmatrix}
\mathbf{P}_{i-p,j-q} & \cdots & \mathbf{P}_{i-p,j}\\
\vdots& \cdots &\vdots\\
\mathbf{P}_{i,j-q} & \cdots & \mathbf{P}_{i,j}\\
\end{pmatrix}
\begin{pmatrix}
N_{j-q,q}(v)\\
\vdots\\
N_{j-q,q}(v)\\
\end{pmatrix}_{(q+1)\times 1}
$$

where,

$$\mathbf P_{i,j}=\{x_{i,j},y_{i,j},z_{i,j}\}$$ or

$$\mathbf P_{i,j}=\{x_{i,j},y_{i,j},z_{i,j}, w_{i,j}\}$$

So the middle matrix is actually a tensor with rank 3, like this:

{{{1,1,0.09986},{1,2,0.07914},{1,3,0.80686},{1,4,0.68846},{1,5,0.297402}},
 {{2,1,0.12262},{2,2,0.41283},{2,3,0.38565},{2,4,0.05670},{2,5,-0.12276}},
 {{3,1,0.08301},{3,2,0.13181},{3,3,0.36565},{3,4,0.74432},{3,5,-0.62065}},
 {{4,1,0.12755},{4,2,0.06099},{4,3,0.74465},{4,4,0.18402},{4,5,0.336987}},
 {{5,1,0.12346},{5,2,0.97057},{5,3,0.72663},{5,4,0.23481},{5,5,0.968757}}}

Here is my algorithm and implementation with C code:

```
void calc_bspline_surf(double P, double U, double V, int p, int q, double u, double v, int const dims, double *S) {
double Nu = (double ) malloc((p + 1)*sizeof(double));
double Nv = (double ) malloc((q + 1)*sizeof(double));
double temp = (double ) malloc(dims[2]*sizeof(double));
int m, n, c;
int i, j;
int k, l, r;
int uidx, vidx, idx;
m = dims[0];
n = dims[1];
c = dims[2];
i = find_span_index(p, U, m - 1, u);
j = find_span_index(q, V, n - 1, v);
b_spline_basis(p, U, i, u, Nu);
b_spline_basis(q, V, j, v, Nv);
/main implementation/
uidx = i - p;
for (l = 0; l <= q; l++){
for (r = 0; r <= c - 1; r++){
temp[r] = 0.0;
}
vidx = j - q + l;
for (k = 0; k <= p; k++){
idx = cn(uidx + k) + c*vidx;
for (r = 0; r <= c - 1; r++){
temp[r] = temp

Readability & Maintainability

-
Right now your code is a bit hard to read. A way to greatly help this would be to use better variable names.

-
Put the variable declarations to separate lines and initialize them to some value at the same time. From Code Complete, 2d Edition, p. 759:


With statements on their own lines, the code reads from top to bottom,
instead of top to bottom and left to right. When you’re looking for a
specific line of code, your eye should be able to follow the left
margin of the code. It shouldn’t have to dip into each and every line
just because a single line might contain two statements.

-
You should declare your for loops as such:

for(int i = 0; i < size; ++i)

Note that this was introduced in the C99 standard. There is no reason you should not be using this standard in your code.

Design

-
Group your variables in structs to make them more maintainable, a
good hint that you should be doing this is when you are passing in
9 parameters to your function.

-
Return the output vector; it makes your function easier to use.
You'll have to allocate the space in your method and free it
elsewhere. This also eliminates the potential bug discussed below.

-
Don't cast malloc().

-
Right now you are using this loop to zero out temp (or at least part of it)

for (r = 0; r <= c - 1; r++){
    temp[r] = 0.0;
}

The original algorithm doesn't implement this in a separate loop, and I don't recommend you do either. It's a loop in the code you don't need.

-
These lines:

temp[r] = temp[r] + Nu[k]*P[idx+r];
S[r] = S[r] + Nv[l]*temp[r];

Can be shortened to:

temp[r] += Nu[k]*P[idx+r];
S[r] += Nv[l]*temp[r];

-
Set those pointers you free()ed to NULL right after. This prevents us from trying to access the data afterwards, which will result in undefined and usually catastrophic behavior.

Bug

-
Right now you are relying on the passed in values already stored in S[r] when you add Nv[l]*temp[r] to them. In the original NURBS A3.5 algorithm, they are 0'ed out right before their use. You could overwrite all these values to 0 with a memset() in our function and continue as is. However, overwriting something the user passes when we could've started from scratch is usually not expected behavior.

Create the array at the top of the function like you do with the other arrays, using calloc().

StackExchange Code Review Q#134614, answer score: 8