c++ - 如何在 Mexfile 中的 matlab(矩阵、单元格)和 c++(向量或自定义类)之间正确转换变量
问题描述
我刚开始使用mex文件,在mexfunction中转换变量时遇到了一些问题,特别是当我想将矩阵转换为自定义类(例如“vema.h”中的Vector)和单元格时MexFunction 中的向量。我们已经提出了这个问题(链接:MATLAB crash when I run MEX file),在这里我们重新组织一下。
下面是计算函数:
#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;
}
功能:
/* 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);
}
}
下面是在 Matlab 中调用 Mexfile 的部分:
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;
以及 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) {};
...
}
谢谢你。
小玉。
解决方案
让我向您展示一个简单的示例,说明如何使用以下方法将数据从 MATLAB 获取到 C++ 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]);
}
但是,请注意,这仅在您的输入mexfucntion
是 afloat
或在 MATLAB 中时才有效single
。
在 MATLAB 中:
A=rand(10,1);
mymexfucntion(A);
A
因为是 type会导致崩溃double
。反而
A=rand(10,1);
mymexfucntion(single(A));
将工作。因为它与 C++ 变量的数据类型相同。如果 C++ 变量是类、向量或简单类型,则无关紧要,如果要强制转换mxGetData
,MATLAB 变量需要是相同的类型。你不能这样做,Vector
因为你编造了。
这一切对您的问题意味着什么?简单:您需要编写自己的 MATLAB-Vector 转换器。
附带说明:mexGetData
不会给您数据的副本。它会给你一个指向数据的指针。如果你修改它,MATLAB中的变量也会被修改,小心!
推荐阅读
- python - 我如何模拟 Python “ModuleNotFoundError”
- arcgis - 使用 JavaScript API 在 ArcGIS Online 上访问安全的 FeatureLayer
- php - Laravel 8 - 未找到基表或视图:1146 表“laravel8.brand”不存在
- go - 如何在 go 中的 map 中获取 for 循环中的所有数据?
- c# - 从 dotnet (Python.NET) 取消 python 函数
- javascript - 使用网格模板自动调整最后一行以占据整列
- flutter - 如何在 SafeArea 中包装 showModalBottomSheet?- 颤振
- c - C 由于预增导致的错误输出
- reactjs - 网络错误 axios 试图从 API 获取数据
- mysql - 如何使这个基本的 MySQL SELECT 查询在 10 亿行上更快?