I've just begun to use mex file, and I met some problems while converting varibales in the mexfunction, especially when I want to convert a matrix to a self-defining class (e.g. Vector in "vema.h") and a cell to a vector in MexFunction. We have posed this question (link:MATLAB crashes when I run MEX file), here we re-organize it.
Below is computational function:
#include "mex.h"
#include "matrix.h"
#include <omp.h>
#include "eig3.h"
#include <stdlib.h>
#include <iostream>
#include <cmath>
#include <vector>
#include <fstream>
#include <iomanip>
#include <sys/types.h>
#include <sys/stat.h>
#include "vema.h" //in vema.h, 'Vector' class has been defined inside which is used in the definition of Fb and V in mexfunction.
using namespace std;
Vector closestPointTriangle(Vector&, Vector&, Vector&, Vector&, double&, double&, double&);
// Generates pomwSize-triangle proximity lists using the linked cell algorithm
void createNNLtriangle(vector<int>* NNLt, Vector* Ut, Vector* faces, int* SN, mwSize nsn, mwSize nf, double hs, double bw, double mw) {
int mx = max(1, (int)(bw/mw)); // ** = 40 cells bw=3.2, mw=0.08
vector<int> head(mx*mx*mx, -1);
vector<int> list(nf);
// std::vector<int> head(mx*mx*mx, -1); //****** mx*mx*mx cells nomber, size mx*mx*mx vector with all values are -1, 40*40*40 = 64000
// std::vector<int> list(nf); // **** nf = 101882
int xa, ya, za, xi, yi, zi;
double ub, vb, wb;
int pt, tri;
Vector cog;
for (int i = 0; i < nf; i++) { // Divide triangle faces mwSizeo cells, i index of face
//Vector cog = (Ut[faces[i].n1] + Ut[faces[i].n2] + Ut[faces[i].n3])/3.0;
cog = (Ut[(int)faces[i].x] + Ut[(int)faces[i].y] + Ut[(int)faces[i].z])/3.0;
int xa = (int)((cog.x + 0.5*bw)/bw*mx);
int ya = (int)((cog.y + 0.5*bw)/bw*mx);
int za = (int)((cog.z + 0.5*bw)/bw*mx);
int tmp = mx*mx*za + mx*ya + xa; // *** 1641838 > 64000
list[i]=head[mx*mx*za + mx*ya + xa];
head[mx*mx*za + mx*ya + xa] = i;
}
#pragma omp parallel for
for (int i = 0; i < nsn; i++) { // Search cells around each pomwSize and build proximity list
int pt = SN[i];
NNLt[i].clear();
int xa = (int)((Ut[pt].x + 0.5*bw)/bw*mx);
int ya = (int)((Ut[pt].y + 0.5*bw)/bw*mx);
int za = (int)((Ut[pt].z + 0.5*bw)/bw*mx);
for (int xi = max(0, xa-1); xi <= min(mx-1, xa+1); xi++)// *** Browse head list
for (int yi = max(0, ya-1); yi <= min(mx-1, ya+1); yi++)
for (int zi = max(0, za-1); zi <= min(mx-1, za+1); zi++) {
int tri = head[mx*mx*zi + mx*yi + xi];
while (tri != -1) {
if ( pt != (int)faces[tri].x && pt != (int)faces[tri].y && pt != (int)faces[tri].z ) {
if ( (closestPointTriangle(Ut[pt], Ut[(int)faces[tri].x], Ut[(int)faces[tri].y], Ut[(int)faces[tri].z], ub, vb, wb) - Ut[pt]).length() < hs) {
NNLt[i].push_back(tri);
}
}
tri = list[tri];
}
}
}
}
// Returns the closest pomwSize of triangle abc to pomwSize p ***** a or b or c, if not, pt projection through the barycenter inside the triangle
Vector closestPointTriangle(Vector& p, Vector& a, Vector& b, Vector& c, double& u, double& v, double& w) {
Vector ab = b - a;
Vector ac = c - a;
Vector ap = p - a;
double d1 = ab.dot(ap);
double d2 = ac.dot(ap);
if (d1 <= 0.0 && d2 <= 0.0) {
u = 1.0;
v = 0.0;
w = 0.0;
return a;
}
Vector bp = p - b;
double d3 = ab.dot(bp);
double d4 = ac.dot(bp);
if (d3 >= 0.0 && d4 <= d3) {
u = 0.0;
v = 1.0;
w = 0.0;
return b;
}
double vc = d1*d4 - d3*d2;
if (vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0) {
v = d1 / (d1 - d3);
u = 1.0 - v;
w = 0.0;
return a + ab * v;
}
Vector cp = p - c;
double d5 = ab.dot(cp);
double d6 = ac.dot(cp);
if (d6 >= 0.0 && d5 <= d6) {
u = 0.0;
v = 0.0;
w = 1.0;
return c;
}
double vb = d5*d2 - d1*d6;
if (vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0) {
w = d2 / (d2 - d6);
u = 1.0 - w;
v = 0.0;
return a + ac * w;
}
double va = d3*d6 - d5*d4;
if (va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0) {
w = (d4 - d3) / ((d4 - d3) + (d5 - d6));
u = 0.0;
v = 1.0 - w;
return b + (c - b) * w;
}
double denom = 1.0 / (va + vb + vc);
v = vb * denom;
w = vc * denom;
u = 1.0 - v - w;
return a + ab * v + ac * w;
}
mexfunction :
/* The gateway function */
void mexFunction(mwSize nlhs, mxArray *plhs[],
mwSize nrhs, const mxArray *prhs[])
{
vector<int> *NNLt;
Vector *V;
Vector *Fb;
int *sn;
mwSize nsn;
mwSize nf;
double hs;
double bw;
double mw;
mwSize i;
/* check for proper number of arguments */
if(nrhs!=9) {
mexErrMsgIdAndTxt("MyToolbox:arrayProduct:nrhs","Nine inputs required.");
}
if(nlhs!=1) {
mexErrMsgIdAndTxt("MyToolbox:arrayProduct:nlhs","One output required.");
}
/* get the value of the scalar input */
nsn = mxGetScalar(prhs[4]);
nf = mxGetScalar(prhs[5]);
hs = mxGetScalar(prhs[6]);
bw = mxGetScalar(prhs[7]);
mw = mxGetScalar(prhs[8]);
/* create a pomwSizeer to the real data in the input matrix */
V = (Vector *)mxGetData(prhs[1]);
Fb = (Vector *)mxGetData(prhs[2]);
sn = (int *)mxGetData(prhs[3]);
for(i=0; i<nsn; i++) {
mxArray *tmp = mxGetCell(prhs[0],i);
std::size_t N = mxGetN(tmp); // number of columns
vector<std::vector<int>> NNLt(nsn, std::vector<int>(N));
double* ptr = mxGetPr(tmp);
copy(ptr, ptr+N, NNLt[i].begin());
};
/* call the computational routine */
createNNLtriangle(NNLt, V, Fb, sn, nsn, nf, hs, bw, mw);
/* create the output matrix */
plhs[0] = mxCreateCellMatrix(nsn,50);
/* get a pomwSizeer to the real data in the output matrix */
for(i=0;i<nsn;i++){
mxArray* tmp = mxCreateDoubleMatrix(1, NNLt[i].size(), mxREAL);
copy(NNLt[i].begin(), NNLt[i].end(), mxGetPr(tmp));
mxSetCell(plhs[0], i, tmp);
}
}
Below is the part to call Mexfile in Matlab:
NNLt = cell(1,nsn);
V = A(1:n,1:3); //A: double Matrix
Fb = A(n:2*n,1:3);
nf = size(Fb,1);
sn = B(1,:); //B: int Matrix
mw = 8.0*a;
hs = 0.6*a;
hc = 0.2*a;
nsn=length(sn);
parfor i = 1:nsn
maxDist = max(maxDist, length(V(sn(i),:) - Vtold(i,:)));
end;
if maxDist > 0.5*(hs-hc)
[NNLt] = createNNLtriangle(NNLt, V, Fb, sn, nsn, nf, hs, bw, mw);
for i = 1:nsn
Vtold(i,:) = V(sn(i),:);
end;
end;
And a part of vema.h :
#include <cmath>
class Vector{
public:
double x, y, z;
Vector(): x(0.0), y(0.0), z(0.0) {};
Vector(double ax, double ay, double az): x(ax), y(ay), z(az) {};
...
}
Thank you.
Xiaoyu.
Let me show you a simple example of how to get the data from MATLAB to C++ using mex:
void mexFunction(int nlhs , mxArray *plhs[],
int nrhs, mxArray const *prhs[])
{
float * myvar= static_cast<float *>(mxGetData(prhs[0]));
mexPrintf("%f \n",myvar[0]);
}
However, note, this will only work if the input to your mexfucntion is a float, or in MATLAB, single.
In MATLAB:
A=rand(10,1);
mymexfucntion(A);
Will cause a crash because A is type double. Instead
A=rand(10,1);
mymexfucntion(single(A));
will work. Because its the same data type as the C++ variable. It is irrelevant if the C++ variable is a class, a vector or a simple type, if you want to cast and mxGetData it, the MATLAB variable needs to be of the same type. You can not do that with Vector because you made that up.
What does all this means to your problem? Simple: You need to write your own MATLAB-Vector converter.
As a side note: mexGetData will not give you a copy of the data. It will give you a pointer to the data. If you modify it, the variable in MATLAB will also be modified, caution!