# Two knot-removal function for curve and surface

Recently, I have been using the Wolfram LibraryLink wrapper. By that technique, I could call the function written in C from Wolfram Mathematica.

Here are two C wrapper functions, BSplineCurveKnotsRemoveAll() and BSplineSurfaceKnotsRemoveAll():

• BSplineCurveKnotsRemoveAll(P,U,p,TOL) removes as many knots as possible of a B-spline curve according to the tolerance TOL and remove_all_curve_knots().

• BSplineSurfaceKnotsRemoveAll(P,U,p,TOL) removes as many knots as possible of a B-spline surface in the row direction via the the tolerance TOL and remove_all_surf_row_knots().

For the curve case, the control points is a vector, here $\mathbf{P}_i=\{x_i, y_i\}$:

$$\{\mathbf{P}_0, \mathbf{P}_1, \ldots, \mathbf{P_n}\}$$

While for surface case, the control nets is a matrix:

$$\begin{Bmatrix} \\ \{\mathbf{P}(0,0), \mathbf{P}(0,1), \space \ldots, \mathbf{P}(0,n)\}, \\ \{\ldots, \space \ldots, \space \ldots, \space \ldots\}, \\ \{\mathbf{P}(m,0), \mathbf{P}(m,1), \space \ldots, \mathbf{P}(m,n)\} \\ \end{Bmatrix}$$

where $\mathbf{P}_{i,j}=\{x_{i,j},y_{i,j},z_{i,j}\}$

BSplineCurveKnotsRemoveAll

//BSplineCurveKnotsRemoveAll(P,U,p,TOL)

DLLEXPORT int BSplineCurveKnotsRemoveAll(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res) {
MTensor tensor_P, tensor_U;
mreal *P, *U;
mint p;
mreal TOL;
mint const *P_dims;
CURVE curve;
mint t;
int err;
mint dims, n, c;
MTensor tensor_Pw;
mreal *Pw;
mint i, j, idx;
//get the arguments P, U, p, TOL
tensor_P = MArgument_getMTensor(Args);
P = libData->MTensor_getRealData(tensor_P);
tensor_U = MArgument_getMTensor(Args);
U = libData->MTensor_getRealData(tensor_U);
p = MArgument_getInteger(Args);
TOL = MArgument_getReal(Args);
P_dims = libData->MTensor_getDimensions(tensor_P);
n = P_dims - 1;
c = P_dims;
//assign the corresponding value to CURVE struct variables
curve.P = P;
curve.U = U;
curve.p = p;
curve.dims = P_dims;
/*call the function to remove as many knots as possible according to
the tolerance TOL, and t is the number of knots that be removed*/
t = remove_all_curve_knots(&curve, TOL);
/* t==0, i.e., no knot is removed, return the control points P directly,
otherwise, copy the valid control points to Pw*/
if (t == 0) {
MArgument_setMTensor(Res, tensor_P);
}
else{
/*build the tensor*/
dims = n - t + 1;
dims = c;
err = libData->MTensor_new(MType_Real, 2, dims, &tensor_Pw);
Pw = libData->MTensor_getRealData(tensor_Pw);
for (i = 0; i < dims; i++) {
for (j = 0; j < c; j++) {
idx = c * i + j;
Pw[idx] = P[idx];
}
}
MArgument_setMTensor(Res, tensor_Pw);
}
return LIBRARY_NO_ERROR;
}


BSplineSurfaceKnotsRemoveAll

//BSplineSurfaceKnotsRemoveAll(P,U,p,TOL)

DLLEXPORT int BSplineSurfaceKnotsRemoveAll(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res) {
MTensor tensor_P, tensor_U;
mreal *P, *U;
mint p;
mreal TOL;
mint const *P_dims;
CURVE curves;
mint t;
int err;
mint dims, m, n, c;
MTensor tensor_Pw;
mreal *Pw;
mint i, j, k, idx;
//get the arguments P, U, p, TOL
tensor_P = MArgument_getMTensor(Args);
P = libData->MTensor_getRealData(tensor_P);
tensor_U = MArgument_getMTensor(Args);
U = libData->MTensor_getRealData(tensor_U);
p = MArgument_getInteger(Args);
TOL = MArgument_getReal(Args);
P_dims = libData->MTensor_getDimensions(tensor_P);
m = P_dims - 1;
n = P_dims - 1;
c = P_dims;
//assign the corresponding value to CURVE struct variables
curves.P = P;
curves.U = U;
curves.p = p;
curves.dims = P_dims;
/*call the function to remove as many knots as possible according to
the tolerance TOL, and t is the number of knots that be removed*/
t = remove_all_surf_row_knots(&curves, TOL);
/* t==0, i.e., no knot is removed, return the control nets P directly,
otherwise, copy the valid control nets to Pw*/
if (t == 0) {
MArgument_setMTensor(Res, tensor_P);
}
else{
/*build the tensor*/
dims = m - t + 1;
dims = n + 1;
dims = c;
err = libData->MTensor_new(MType_Real, 3, dims, &tensor_Pw);
Pw = libData->MTensor_getRealData(tensor_Pw);
for (i = 0; i < dims; i++) {
for (j = 0; j < dims; j++) {
for (k = 0; k < c; k++) {
idx = c * (n + 1) * i + c * j + k;
Pw[idx] = P[idx];
}
}
}
MArgument_setMTensor(Res, tensor_Pw);
}
return LIBRARY_NO_ERROR;
}


In addition, the definition of struct CURVE:

typedef struct curv{
double *P;              //control points of curve or curve group
double *U;              //knot vector
int p;                  //degree of curve
const size_t *dims;     //dimentions of curve or curve group
} CURVE;


I think the two C functions are similar in some sense, apart from the rank of the argument P. For the curve case, the argument P is a tensor with rank $2$ and dimensions {n + 1, c}. While for a surface case, the variable P is a tensor with rank $3$ and dimensions {m + 1, n + 1, c}.

Is it possible to use a C function BSplineKnotsRemoveAll() to merge them?

It's difficult to address the question you're most interested in without more comments (or documentation for the library you're using). However, the differences between the two functions do seem very minor (especially once I rename curves to curve in the second one).

There are three comments in each function:

    //implementation start
/*main implementation*/
/*build the tensor*/


The third one is mildly useful, but the other two aren't really. What would be useful is some indication of the meanings of the 18 or so variables.

In particular, you mention that the dimensions of the tensors are {n + 1, c} and {m + 1, n + 1, c}, but there's no obvious reason in the code for using n + 1 instead of n, or m + 1 instead of m. Both are assigned as a dimension of an input tensor minus one, and then are only ever used in expressions in which they're incremented. I find this quite confusing.

### Differences

There are fundamentally two differences between the functions. The first is which function they call to remove knots, which would obviously be a function pointer parameter of the function you propose to extract. The second is the /*build the tensor*/ section. The latter is not as hard to commonalise as it seems on first sight. Supposing that you have a rank variable (either a parameter or obtained via something like - I'm guessing - libData->MTensor_getRank(tensor_P)). Then the common code would seem to be

        /*build the tensor*/
memcpy(&dims, P_dims, rank * sizeof(mint));
dims -= t;
err = libData->MTensor_new(MType_Real, rank, dims, &tensor_Pw);
Pw = libData->MTensor_getRealData(tensor_Pw);
dimProd = dims;
for (i = 1; i < rank; i++) {
dimProd *= dims[i];
}
memcpy(Pw, P, dimProd * sizeof(mreal));
MArgument_setMTensor(Res, tensor_Pw);


I'm making a few assumptions here. In particular, I'm assuming that it's necessary to copy P_dims because the original's values must not be destroyed.

### Error handling

I've left the error handling in the refactored code above as it was in the original, but it's really bad form to ignore the error code from libData->MTensor_new. If the allocation fails, you're going to try to access an invalid memory address.

• About the reason in the code for using n + 1 instead of n,please see my Update. In addition, you are right, there indeed is a C function called libData->MTensor_getRank() – xyz Aug 24 '16 at 1:23
• For the variable dims in code memcpy(&dims, P_dims, rank * sizeof(mint));, I don't know how to declare it. mint dims? or mint dims? or mint *dims? – xyz Aug 24 '16 at 10:50
• This is one of the things I'm not clear about with respect to the library: is it actually necessary to copy in the first place, or can you just modify P_dims? Since the knot removal functions seem to modify curve.P, it's even possible that my suggested common code is buggy because t has already been removed from P_dims. – Peter Taylor Aug 24 '16 at 10:52
• Yes, when calling the knots removal function, curve.P will be changed. While for its dimensions P_dims, which is constant in whole process. – xyz Aug 24 '16 at 11:00
• In addition, I don't know why use &dims, rather than dims? – xyz Aug 24 '16 at 11:08