#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include "NURBS.h"
#include "mex.h"
#define free2Darray
#define dersBasisFuns
#define BasisFuns
#define FindSpan
#define init2DArray
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
/* Return the NURBS first derivatives w.r.t xi and eta to matlab
// All non-zero derivatives at point [xi,eta] are computed.
//
// We expect the function to be called as
// [dRdxi dRdeta] = NURBSinterpolation(xi, p, q, knotU, ...
// knotV, weights)
//
// xi = point, [xi eta], where we want to interpolate
// knotU, knotV = knot vectors
// weights = vector of weights
//
// Vinh Phu Nguyen, nvinhphu@gmail.com
*/
double *xi;
int p_in;
int q_in;
int p;
int q;
double *knotU;
double *knotV;
int numKnotU;
int numKnotV;
int nU;
int nV;
int noFuncs;
double *weight;
int numWeights;
double tol;
double *N;
double *M;
double **dersN;
double **dersM;
int spanU;
int spanV;
int i, j, k, c;
int uind;
int vind;
double w;
double dwdxi;
double dwdet;
double wgt;
double *dRdxi;
double *dRdet;
double fac;
int g;
/* Examine input (right-hand-side) arguments. */
mexPrintf("\nThere are %d right-hand-side argument(s).", nrhs);
for (g=0; g<nrhs; g++) {
mexPrintf("\n\tInput Arg %i is of type:\t%s\n ",g,mxGetClassName(prhs[g]));
}
/* Examine output (left-hand-side) arguments. */
mexPrintf("\n\nThere are %d left-hand-side argument(s).\n", nlhs);
if (nlhs > nrhs)
mexErrMsgIdAndTxt( "MATLAB:mexfunction:inputOutputMismatch",
"Cannot specify more outputs than inputs.\n");
if(nrhs != 6) mexErrMsgTxt("You fool! You haven't passed in 6 arguments to the function."
"We expect it to be in the form [dRdxi dRdeta] = NURBSinterpolation(xi, p, q, knotU, knotV, weights)\n");
/*First get the inputs */
if( !mxIsSingle(prhs[0]) ) {
mexErrMsgTxt("First input must be single");
}
xi = mxGetData(prhs[2]);
p_in = mxGetScalar(prhs[0]);
q_in = mxGetScalar(prhs[1]);
printf("xi1= %.5f\n", xi[0]);
printf("xi2= %.5f\n", xi[1]);
p = (int) p_in;
q = (int) q_in;
printf("p= %d\n", p);
printf("q= %d\n", q);
knotU = mxGetPr(prhs[3]);
knotV = mxGetPr(prhs[4]);
printf("knotU1= %.2f\n", knotU[0]);
printf("knotU2= %.2f\n", knotU[1]);
printf("knotU3= %.2f\n", knotU[2]);
printf("knotU4= %.2f\n", knotU[3]);
printf("knotU5= %.2f\n", knotU[4]);
printf("knotU6= %.2f\n", knotU[5]);
printf("knotU7= %.2f\n", knotU[6]);
numKnotU = mxGetNumberOfElements(prhs[3]);
numKnotV = mxGetN(prhs[4]);
printf("numKnotU= %d\n", numKnotU);
nU = (int)numKnotU - 1 - p - 1;
nV = (int)numKnotV - 1 - q - 1;
printf("nU= %d\n", nU);
printf("nV= %d\n", nV);
noFuncs = (p+1)*(q+1);
printf("noFuncs= %d\n", noFuncs);
weight = mxGetPr(prhs[5]);
numWeights = mxGetN(prhs[5]);
printf("weight2= %.5f\n", weight[1]);
printf("numWeights= %d\n", numWeights);
tol = 100*DBL_EPSILON;
printf("tol= %d\n", tol);
/* if(fabs(xi[0]-knotU[numKnotU-1]) < tol)
xi[0] = knotU[numKnotU-1] - tol;
if(fabs(xi[1]-knotV[numKnotV-1]) < tol)
xi[1] = knotV[numKnotV-1] - tol;
*/
/* and evaluate the non-zero univariate B-spline basis functions
* and first derivatives
*/
N = (p+1);
printf("N= %d\n", N);
M = (q+1);
printf("N= %d\n", M);
dersN = init2DArray((int)nU+1, (int)p+1);
dersM = init2DArray((int)nV+1, (int)q+1);
printf("dersN= %d\n", dersN);
printf("dersM= %d\n", dersM);
spanU = mxGetScalar(FindSpan(3, 2, 0.05640, knotU));
spanV = mxGetScalar(FindSpan(nV, q, xi[1], knotV));
printf("spanU= %d\n", spanU);
printf("spanV= %d\n", spanV);
/* BasisFunsU= BasisFuns ((int)spanU, (int)xi[0], (int)p, (double *)knotU, (double *)N);
BasisFunsV= BasisFuns ((int)spanV,(int) xi[1], (int)q, (double *)knotV, (double *)M);
printf("BasisFunsU= %d\n", BasisFunsU);
printf("BasisFunsV= %d\n", BasisFunsV);
dersBasisFunsU= dersBasisFuns ((int)spanU, (int)xi[0], (int)p, (int)nU, (double *)knotU, (double **)dersN);
dersBasisFunsV= dersBasisFuns ((int)spanV, (int)xi[1], (int)q, (int)nV, (double *)knotV, (double **)dersM);
printf("dersBasisFunsU= %d\n", dersBasisFunsU);
printf("dersBasisFunsV= %d\n", dersBasisFunsV);
/* and create NURBS approximation */
/*
for(i=0;i<=p;i++){
printf("dNdxi= %f\n", dersN[1][i]);
printf("dNdet= %f\n", dersM[1][i]);
}*/
/*
uind = (int)spanU - (int)p;
vind;
w = 0.0;
dwdxi = 0.0;
dwdet = 0.0;
wgt;
for(j = 0; j <= q; j++)
{
vind = spanV - q + j;
for(i = 0; i <= p; i++)
{
c = (int)uind + (int)i + (int)vind * ((int)nU+1);
wgt = weight[c];
w += N[i] * M[j] * wgt;
dwdxi += dersN[1][i] * M[j] * wgt;
dwdet += dersM[1][j] * N[i] * wgt;
}
}
/* create output */
/*
plhs[0] = mxCreateDoubleMatrix(1,noFuncs,mxREAL);
plhs[1] = mxCreateDoubleMatrix(1,noFuncs,mxREAL);
dRdxi = mxGetPr(plhs[0]);
dRdet = mxGetPr(plhs[1]);
uind = spanU - p;
k = 0;
for(j = 0; j <= q; j++)
{
vind = spanV - q + j;
for(i = 0; i <= p; i++)
{
c = uind + i + vind*(nU+1);
fac = weight[c]/(w*w);
dRdxi[k] = (dersN[1][i]*M[j]*w - N[i]*M[j]*dwdxi) * fac;
dRdet[k] = (dersM[1][j]*N[i]*w - N[i]*M[j]*dwdet) * fac;
k += 1;
}
}
/*mexPrintf("\nWe have a knot vector with %d components\n", numWeights);
for(k = 0; k < numWeights; k++) mexPrintf("%2.2f\t", weight[k]);
mexPrintf("\n");*/
/*
free(N);
free(M);
free2Darray(dersN, (nU+1));
free2Darray(dersM, (nV+1));
*/
}
