transform.c

三维张量在不同坐标系统的变换

/*
 * FUNCTION
 * otr = transform(itr,tmx)
 *
 * DESCRIPTION
 * transform 3D-tensor (Euclidean or Cartesion tensor) of any order (>0) to another coordinate system
 *
 * PARAMETERS
 * otr = output tensor, after transformation; has the same dimensions as the input tensor
 * itr = input tensor, before transformation; should be a 3-element vector, a 3x3 matrix, or a 3x3x3x... multidimensional array, each dimension containing 3 elements
 * tmx = transformation matrix, 3x3 matrix that contains the direction cosines between the old and the new coordinate system
 */

#include "mex.h"

void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )
     
{ 
  int id, ie, oe;
  int nd, ne;
  int *iie, *ioe, *cne;
  double *itr, *tmx, *otr;
  double pmx;
    
  itr = mxGetPr(prhs[0]);                       /* get the input tensor */
  tmx = mxGetPr(prhs[1]);                       /* get the transformation matrix */
  
  ne = mxGetNumberOfElements(prhs[0]);          /* number of tensor elements */
  nd = mxGetNumberOfDimensions(prhs[0]);        /* number of tensor dimensions, i.e. order of tensor */
  nd = ((ne==3) ? 1 : nd);                      /* order of tensor is 1 in case of a 3x1 or 1x3 vector */
 
  plhs[0] = mxCreateNumericArray(nd, mxGetDimensions(prhs[0]), mxDOUBLE_CLASS, mxREAL);   /* create and initialise the output array */
  otr = mxGetPr(plhs[0]);                       /* set the output tensor */

  iie = mxCalloc(nd, sizeof(int));              /* initialise vector with indices of input tensor element */
  ioe = mxCalloc(nd, sizeof(int));              /* initialise vector with indices of output tensor element */
  cne = mxCalloc(nd, sizeof(int));              /* initialise vector with cumulative number of elements for each dimension (divided by three) */
  cne[0] = 1;                                   /* create the vector */
  for (id=1; id<nd; id++) {
     cne[id] = cne[id-1] * 3;
  }
  
  for (oe=0; oe<ne; oe++) {                     /* loop over all output elements */
     for (id=0; id<nd; id++) {
        ioe[id] = ((oe/cne[id]) % 3);           /* calculate indices of current output tensor element */
     }
     for (ie=0; ie<ne; ie++) {                  /* loop over all output elements */
        pmx = 1;                                /* initialise product of transformation matrices */
        for (id=0; id<nd; id++) {               /* loop over all dimensions */
           iie[id] = ((ie/cne[id]) % 3);        /* calculate indices of current input tensor element */
           pmx = pmx * tmx[ioe[id]+3*iie[id]];  /* create product of transformation matrices */
        }
        otr[oe] = otr[oe] + pmx * itr[ie];      /* add product of transformation matrices and input tensor element to output tensor element */
     }
  }

}