Hi Guys!

Can eigen work with 3d matrices? Im currently rewriting matlab code in c++ and i need to implement the meshgrid(x, y, z) function.

Can anyone help me out? Would be really appreciated

Eigen does not support multi-dimensional arrays or tensors. So you'll have to mimic it with a vector of 2D matrices, and setLinSpaced.

I'll post the matlab code and maybe you'll see the problem im having

Code: Select all

nn=-n:1:n;                            % Index for the sequence
rms=nn+0.5-0.5*(-1).^nn;              % Part of equations 2,3,& 4
srcs=(-1).^(nn);                      % part of equations 2,3,& 4
xi=srcs*src(1)+rms*rm(1)-mic(1);      % Equation 2
yj=srcs*src(2)+rms*rm(2)-mic(2);      % Equation 3
zk=srcs*src(3)+rms*rm(3)-mic(3);      % Equation 4

[i,j,k]=meshgrid(xi,yj,zk);           % convert vectors to 3D matrices
d=sqrt(i.^2+j.^2+k.^2);               % Equation 5
time=round(fs*d/343)+1;               % Similar to Equation 6
             
[e,f,g]=meshgrid(nn, nn, nn);         % convert vectors to 3D matrices
c=r.^(abs(e)+abs(f)+abs(g));          % Equation 9
e=c./d;                               % Equivalent to Equation 10

h=full(sparse(time(:),1,e(:)));       % Equivalent to equation 11
h=h/max(abs(h));                      % Scale output



I dont know how i would calculate the euclidean norm ( d=sqrt(i.^2+j.^2+k.^2) ) in this way.... Seems like working with a vector of matrices wouldnt be very friendly. Also the code needs to run as fast as possible, the ideal being almost in realtime.

Not tested, but that's the general idea:

Code: Select all

ArrayXf x, y, z;
ArrayXf x2y2(N*N), d(N*N*N);
ArrayXXf::Map(x2y2.data(), N, N) = x.square().replicate(1,N) + y.square().transpose().replicate(N,1);
ArrayXXf::Map(d.data(), N*N, N) = x2y2.replicate(1,N) + z.square().transpose().replicate(N,1);


Thanks i'll give that a go tomorrow!

I've been looking at your code and i dont think this does the same thing as matlabs meshgrid function. I'm really lost

Edit:

I understand that i want to create an array of matrices for this, but i dont understand your code. what is x2y2? And why is d a vector sized N^3, when its supposed to be a Matrix sized N^3? I kind of imagine working with a variable like this:

NxN
NxN
NxN
v = ( . )
.
.
NxN

So v is an Array of matrices. Is this possible in eigen?

Since you only apply coefficient-wise operations on d, c, r, etc. it is finally simpler to represent them as 1D arrays. Moreover, at this end your 3D matrices are linearized anyway.

x2y2 is a N^2 array where x2y2(i+j*N) = x_i^2 + y_j^2.

The same technique is used to compute d from x2y2. And you keep using the same approach to compute c. The other operations are straightforward.

Ok I understand now and could extrapolate from this, but im hitting a breakpoint when i debug the application for some reason and i dont know why...

Heres the code, i know it could be more efficient, but as i said im pretty new to all this and am trying my best

Code: Select all

#include <iostream>
#include <Eigen/Dense>

using namespace Eigen;
using namespace std;

int main()
{

   float fSamplerate = 44100;

   double Xmic = 1;
   double Ymic = 1;
   double Zmic = 1;

   double Xsrc = 3;
   double Ysrc = 2;
   double Zsrc = 1;

   double Xrm = 5;
   double Yrm = 4;
   double Zrm = 3;

   int n = 5;
   int ntemp,icounter;
   double dTemp;
   double r = 0.4;

   ArrayXd nn(2 * n + 1);

   int N = nn.rows();

   ArrayXd rms(N);
   ArrayXd srcs(N);

   ArrayXd xi(N);
   ArrayXd yj(N);
   ArrayXd zk(N);

   ArrayXd x2y2(N*N);
   ArrayXd d(N*N*N);


   // init nn
   icounter = 0;
   for(ntemp = -n; ntemp <= n; ntemp++)
   {
      nn(icounter) = ntemp;
      icounter++;
   }

   //init rms and srcs
   for(icounter = 0; icounter <= 2*n; icounter++)
   {
      dTemp = pow(-1.0,nn(icounter));
      rms(icounter) = nn(icounter) + (0.5 - 0.5 * dTemp); 
      srcs(icounter) = dTemp;
   }

   // calculate relative positions along x,y,z axis
   xi = srcs * Xsrc + (rms * Xrm) - Xmic;
   yj= srcs * Ysrc + (rms * Yrm) - Ymic;
   zk= srcs * Zsrc + (rms * Zrm) - Zmic;

   // calculate d
   ArrayXXd::Map(x2y2.data(), N, N) = xi.square().replicate(1,N) + yj.square().transpose().replicate(N,1);
   ArrayXXd::Map(d.data(), N*N, N) = x2y2.replicate(1,N) + zk.square().transpose().replicate(N,1);


   cout << "d = " << d << endl;


   std::getchar();
}


I would be really thankful for some help!

ArrayXXf::Map(x2y2.data(), N, N) = x.square().replicate(1,N) + y.square().transpose().replicate(N,1);

when you replicate x and y they become matrices and then you map them to a vector again. How does this work? I'm pretty sure the error im getting is to do with vector dimensions not matching up!

Thank you so much!

It's "zk.square().transpose().replicate(N*N,1)" and not "zk.square().transpose().replicate(N,1)".

Yep thanks! I actually managed to work it out yesterday too, but it took quite a while. So for future reference:

Code: Select all

// Matlab code
// [i,j,k]=meshgrid(xi,yj,zk);
// d=sqrt(i.^2+j.^2+k.^2);

// C++ code

ArrayXd x2y2(N*N);
ArrayXd d(N*N*N);

ArrayXXd::Map(x2y2.data(), N, N) = xi.square().transpose().replicate(N,1) + yj.square().replicate(1,N);
ArrayXXd::Map(d.data(), N*N, N) = x2y2.replicate(1,N) + zk.square().transpose().replicate(N*N,1);

d = d.sqrt();



Thanks again for the help!

btw, if at the end you can come up with some speed comparison with the initial MatLab code, that would be nice!

I'll see if i get around to it and will post results if i do!