Xmipp  v3.23.11-Nereus
Classes | Functions
Multidimensional Arrays
Collaboration diagram for Multidimensional Arrays:

Classes

class  MultidimArray< T >
 
class  MultidimArrayGeneric
 

Functions

template<typename T >
void coreArrayByScalar (const MultidimArray< T > &op1, const T &op2, MultidimArray< T > &result, char operation)
 
template<typename T >
void coreScalarByArray (const T &op1, const MultidimArray< T > &op2, MultidimArray< T > &result, char operation)
 
template<typename T >
void coreArrayByArray (const MultidimArray< T > &op1, const MultidimArray< T > &op2, MultidimArray< T > &result, char operation)
 
template<typename T >
void selfCoreArrayByArrayMask (const MultidimArray< T > &op1, const MultidimArray< T > &op2, MultidimArray< T > &result, char operation, const MultidimArray< T > *mask)
 
template<>
std::ostream & operator<< (std::ostream &ostrm, const MultidimArray< std::complex< double > > &v)
 
template<>
bool operator== (const MultidimArray< std::complex< double > > &op1, const MultidimArray< std::complex< double > > &op2)
 

Functions for all multidimensional arrays

template<typename T1 , typename T2 >
void typeCast (const MultidimArray< T1 > &v1, MultidimArray< T2 > &v2)
 
template<typename T1 >
void typeCastComplex (const MultidimArray< T1 > &v1, MultidimArray< std::complex< double > > &v2)
 
template<typename T1 >
void typeCast (const MultidimArray< T1 > &v1, MultidimArray< T1 > &v2)
 
template<typename T >
void typeCast (const MultidimArray< T > &v1, Matrix1D< T > &v2)
 
template<typename T >
bool operator== (const MultidimArray< T > &op1, const MultidimArray< T > &op2)
 
template<typename T >
bool operator!= (const MultidimArray< T > &op1, const MultidimArray< T > &op2)
 
template<typename T >
void cutToCommonSize (MultidimArray< T > &V1, MultidimArray< T > &V2)
 
void sincos (const MultidimArray< double > &x, MultidimArray< double > &s, MultidimArray< double > &c)
 
void planeFit (const MultidimArray< double > &z, const MultidimArray< double > &x, const MultidimArray< double > &y, double &p0, double &p1, double &p2)
 
template<typename T >
void mod (const MultidimArray< T > &x, MultidimArray< T > &m, double y)
 
template<typename T >
std::ostream & operator<< (std::ostream &ostrm, const MultidimArray< T > &v)
 
template<typename T >
void window2D (const MultidimArray< T > &Ibig, MultidimArray< T > &Ismall, size_t y0, size_t x0, size_t yF, size_t xF)
 
template<typename T >
double correlationIndex (const MultidimArray< T > &x, const MultidimArray< T > &y, const MultidimArray< int > *mask=NULL, MultidimArray< double > *Contributions=NULL)
 

MultidimArrayGenericSpeedUp Speed up macros

#define MULTIDIM_ARRAY_BASE(v)   (*((v).data->im))
 
#define MULTIDIM_ARRAY_GENERIC(v)   (*((v).data))
 
#define MULTIDIM_ARRAY_TYPE(v, type)   (*((MultidimArray<type>*)(v.im)))
 

Detailed Description

Macro Definition Documentation

◆ MULTIDIM_ARRAY_BASE

#define MULTIDIM_ARRAY_BASE (   v)    (*((v).data->im))

Array access.

This macros gives you access to the array (T **)

Definition at line 97 of file multidim_array_generic.h.

◆ MULTIDIM_ARRAY_GENERIC

#define MULTIDIM_ARRAY_GENERIC (   v)    (*((v).data))

Definition at line 101 of file multidim_array_generic.h.

◆ MULTIDIM_ARRAY_TYPE

#define MULTIDIM_ARRAY_TYPE (   v,
  type 
)    (*((MultidimArray<type>*)(v.im)))

Definition at line 105 of file multidim_array_generic.h.

Function Documentation

◆ coreArrayByArray()

template<typename T >
void coreArrayByArray ( const MultidimArray< T > &  op1,
const MultidimArray< T > &  op2,
MultidimArray< T > &  result,
char  operation 
)

Core array by array operation.

It assumes that the result is already resized.

Definition at line 2163 of file multidim_array.h.

2166  {
2167  T* ptrResult=NULL;
2168  T* ptrOp1=NULL;
2169  T* ptrOp2=NULL;
2170  size_t n;
2171  // Loop unrolling
2172  const size_t unroll=4;
2173  size_t nmax=unroll*(op1.nzyxdim/unroll);
2174  switch (operation)
2175  {
2176  case '+':
2177  for (n=0, ptrResult=result.data, ptrOp1=op1.data,ptrOp2=op2.data;
2178  n<nmax; n+=unroll, ptrResult+=unroll, ptrOp1+=unroll, ptrOp2+=unroll)
2179  {
2180  *ptrResult = *ptrOp1 + *ptrOp2;
2181  *(ptrResult+1) = *(ptrOp1+1) + *(ptrOp2+1);
2182  *(ptrResult+2) = *(ptrOp1+2) + *(ptrOp2+2);
2183  *(ptrResult+3) = *(ptrOp1+3) + *(ptrOp2+3);
2184  }
2185  for (n=nmax, ptrResult=result.data+nmax, ptrOp1=op1.data+nmax, ptrOp2=op2.data+nmax;
2186  n<op1.zyxdim; ++n, ++ptrResult, ++ptrOp1, ++ptrOp2)
2187  *ptrResult = *ptrOp1 + *ptrOp2;
2188  break;
2189  case '-':
2190  for (n=0, ptrResult=result.data, ptrOp1=op1.data,ptrOp2=op2.data;
2191  n<nmax; n+=unroll, ptrResult+=unroll, ptrOp1+=unroll, ptrOp2+=unroll)
2192  {
2193  *ptrResult = *ptrOp1 - *ptrOp2;
2194  *(ptrResult+1) = *(ptrOp1+1) - *(ptrOp2+1);
2195  *(ptrResult+2) = *(ptrOp1+2) - *(ptrOp2+2);
2196  *(ptrResult+3) = *(ptrOp1+3) - *(ptrOp2+3);
2197  }
2198  for (n=nmax, ptrResult=result.data+nmax, ptrOp1=op1.data+nmax, ptrOp2=op2.data+nmax;
2199  n<op1.zyxdim; ++n, ++ptrResult, ++ptrOp1, ++ptrOp2)
2200  *ptrResult = *ptrOp1 - *ptrOp2;
2201  break;
2202  case '*':
2203  for (n=0, ptrResult=result.data, ptrOp1=op1.data,ptrOp2=op2.data;
2204  n<nmax; n+=unroll, ptrResult+=unroll, ptrOp1+=unroll, ptrOp2+=unroll)
2205  {
2206  *ptrResult = *ptrOp1 * *ptrOp2;
2207  *(ptrResult+1) = *(ptrOp1+1) * *(ptrOp2+1);
2208  *(ptrResult+2) = *(ptrOp1+2) * *(ptrOp2+2);
2209  *(ptrResult+3) = *(ptrOp1+3) * *(ptrOp2+3);
2210  }
2211  for (n=nmax, ptrResult=result.data+nmax, ptrOp1=op1.data+nmax, ptrOp2=op2.data+nmax;
2212  n<op1.zyxdim; ++n, ++ptrResult, ++ptrOp1, ++ptrOp2)
2213  *ptrResult = *ptrOp1 * *ptrOp2;
2214  break;
2215  case '/':
2216  for (n=0, ptrResult=result.data, ptrOp1=op1.data,ptrOp2=op2.data;
2217  n<nmax; n+=unroll, ptrResult+=unroll, ptrOp1+=unroll, ptrOp2+=unroll)
2218  {
2219  *ptrResult = *ptrOp1 / *ptrOp2;
2220  *(ptrResult+1) = *(ptrOp1+1) / *(ptrOp2+1);
2221  *(ptrResult+2) = *(ptrOp1+2) / *(ptrOp2+2);
2222  *(ptrResult+3) = *(ptrOp1+3) / *(ptrOp2+3);
2223  }
2224  for (n=nmax, ptrResult=result.data+nmax, ptrOp1=op1.data+nmax, ptrOp2=op2.data+nmax;
2225  n<op1.zyxdim; ++n, ++ptrResult, ++ptrOp1, ++ptrOp2)
2226  *ptrResult = *ptrOp1 / *ptrOp2;
2227  break;
2228  }
2229  }
int * nmax
int * n

◆ coreArrayByScalar()

template<typename T >
void coreArrayByScalar ( const MultidimArray< T > &  op1,
const T &  op2,
MultidimArray< T > &  result,
char  operation 
)

Core array by scalar operation.

It assumes that the result is already resized.

This function is not ported to Python.

Definition at line 2448 of file multidim_array.h.

2452  {
2453  T* ptrResult=NULL;
2454  T* ptrOp1=NULL;
2455  size_t n;
2456  switch (operation)
2457  {
2458  case '+':
2459  for (n=0, ptrResult=result.data, ptrOp1=op1.data;
2460  n<op1.zyxdim; ++n, ++ptrResult, ++ptrOp1)
2461  *ptrResult = *ptrOp1 + op2;
2462  break;
2463  case '-':
2464  for (n=0, ptrResult=result.data, ptrOp1=op1.data;
2465  n<op1.zyxdim; ++n, ++ptrResult, ++ptrOp1)
2466  *ptrResult = *ptrOp1 - op2;
2467  break;
2468  case '*':
2469  for (n=0, ptrResult=result.data, ptrOp1=op1.data;
2470  n<op1.zyxdim; ++n, ++ptrResult, ++ptrOp1)
2471  *ptrResult = *ptrOp1 * op2;
2472  break;
2473  case '/':
2474  for (n=0, ptrResult=result.data, ptrOp1=op1.data;
2475  n<op1.zyxdim; ++n, ++ptrResult, ++ptrOp1)
2476  *ptrResult = *ptrOp1 / op2;
2477  break;
2478  case '=':
2479  for (n=0, ptrResult=result.data, ptrOp1=op1.data;
2480  n<op1.zyxdim; ++n, ++ptrResult, ++ptrOp1)
2481  *ptrResult = *ptrOp1 == op2;
2482  break;
2483  }
2484  }
int * n

◆ coreScalarByArray()

template<typename T >
void coreScalarByArray ( const T &  op1,
const MultidimArray< T > &  op2,
MultidimArray< T > &  result,
char  operation 
)

Core array by scalar operation.

It assumes that the result is already resized.

This function is not ported to Python.

Definition at line 2608 of file multidim_array.h.

2612  {
2613  T* ptrResult=NULL;
2614  T* ptrOp2=NULL;
2615  size_t n;
2616  for (n=0, ptrResult=result.data, ptrOp2=op2.data;
2617  n<op2.zyxdim; ++n, ++ptrResult, ++ptrOp2)
2618  switch (operation)
2619  {
2620  case '+':
2621  *ptrResult = op1 + *ptrOp2;
2622  break;
2623  case '-':
2624  *ptrResult = op1 - *ptrOp2;
2625  break;
2626  case '*':
2627  *ptrResult = op1 * *ptrOp2;
2628  break;
2629  case '/':
2630  *ptrResult = op1 / *ptrOp2;
2631  break;
2632  }
2633  }
int * n

◆ correlationIndex()

template<typename T >
double correlationIndex ( const MultidimArray< T > &  x,
const MultidimArray< T > &  y,
const MultidimArray< int > *  mask = NULL,
MultidimArray< double > *  Contributions = NULL 
)

correlationIndex nD

Definition at line 3974 of file multidim_array.h.

3978 {
3980 
3981  double retval = 0, aux;
3982  double mean_x, mean_y;
3983  double stddev_x, stddev_y;
3984 
3985  long N = 0;
3986 
3987  if (mask == NULL)
3988  {
3989  x.computeAvgStdev(mean_x, stddev_x);
3990  y.computeAvgStdev(mean_y, stddev_y);
3991  }
3992  else
3993  {
3994  x.computeAvgStdev_within_binary_mask(*mask, mean_x,stddev_x);
3995  y.computeAvgStdev_within_binary_mask(*mask, mean_y,stddev_y);
3996  }
3997  if (ABS(stddev_x)<XMIPP_EQUAL_ACCURACY ||
3998  ABS(stddev_y)<XMIPP_EQUAL_ACCURACY)
3999  return 0;
4000 
4001  // If contributions are desired. Please, be careful interpreting individual
4002  // contributions to the covariance! One pixel value afect others.
4003  if (Contributions != NULL)
4004  {
4006  {
4007  if (mask != NULL)
4008  if (!A3D_ELEM(*mask,k, i, j))
4009  continue;
4010 
4011  aux = (A3D_ELEM(x, k, i, j) - mean_x) * (A3D_ELEM(y, k, i, j) -
4012  mean_y);
4013  A3D_ELEM(*Contributions, k, i, j) = aux;
4014  retval += aux;
4015  ++N;
4016  }
4017 
4018  FOR_ALL_ELEMENTS_IN_ARRAY3D(*Contributions)
4019  A3D_ELEM(*Contributions, k, i, j) /= ((stddev_x * stddev_y) * N);
4020  }
4021  else
4022  {
4023  if (mask==NULL && x.sameShape(y))
4024  {
4026  retval += DIRECT_MULTIDIM_ELEM(x, n)*DIRECT_MULTIDIM_ELEM(y, n);
4027  N=MULTIDIM_SIZE(x);
4028  retval-=N*mean_x*mean_y;
4029  }
4030  else
4031  {
4033  {
4034  if (mask != NULL)
4035  if (!A3D_ELEM(*mask,k, i, j))
4036  continue;
4037 
4038  retval += (A3D_ELEM(x, k, i, j) - mean_x) *
4039  (A3D_ELEM(y, k, i, j) - mean_y);
4040  ++N;
4041  }
4042  }
4043  }
4044 
4045  if (N != 0)
4046  return retval / ((stddev_x * stddev_y) * N);
4047  else
4048  return 0;
4049 }
#define MULTIDIM_SIZE(v)
void computeAvgStdev(U &avg, U &stddev) const
void computeAvgStdev_within_binary_mask(const MultidimArray< int > &mask, double &avg, double &stddev) const
#define i
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
#define A3D_ELEM(V, k, i, j)
#define FOR_ALL_ELEMENTS_IN_ARRAY3D(V)
#define XMIPP_EQUAL_ACCURACY
Definition: xmipp_macros.h:119
bool sameShape(const MultidimArrayBase &op) const
#define FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(v)
#define ABS(x)
Definition: xmipp_macros.h:142
#define DIRECT_MULTIDIM_ELEM(v, n)
#define j
#define SPEED_UP_tempsInt
Definition: xmipp_macros.h:408
#define FOR_ALL_ELEMENTS_IN_COMMON_IN_ARRAY3D(V1, V2)
int * n

◆ cutToCommonSize()

template<typename T >
void cutToCommonSize ( MultidimArray< T > &  V1,
MultidimArray< T > &  V2 
)

Reduce both volumes to a common size.

Search the range of logical indexes for which both volumes have got valid values, and cut both to that size, the corresponding origin is automatically computed.

V1.startingX() = -2;
V1.startingY() = -2;
V1.startingZ() = -2;
V2.startingX() = 0;
V2.startingY() = 0;
V2.startingZ() = 0;
// V1 and V2 range from (0,0,0)=(z,y,x) to (1,1,0)

Definition at line 3866 of file multidim_array.h.

3867 {
3868  int z0 = XMIPP_MAX(STARTINGZ(V1), STARTINGZ(V2));
3869  int zF = XMIPP_MIN(FINISHINGZ(V1), FINISHINGZ(V2));
3870  int y0 = XMIPP_MAX(STARTINGY(V1), STARTINGY(V2));
3871  int yF = XMIPP_MIN(FINISHINGY(V1), FINISHINGY(V2));
3872  int x0 = XMIPP_MAX(STARTINGX(V1), STARTINGX(V2));
3873  int xF = XMIPP_MIN(FINISHINGX(V1), FINISHINGX(V2));
3874 
3875  V1.window(z0, y0, x0, zF, yF, xF);
3876  V2.window(z0, y0, x0, zF, yF, xF);
3877 }
#define XMIPP_MAX(x, y)
Definition: xmipp_macros.h:193
#define FINISHINGX(v)
#define zF
#define z0
#define FINISHINGZ(v)
#define STARTINGX(v)
#define STARTINGY(v)
void window(MultidimArray< T1 > &result, int n0, int z0, int y0, int x0, int nF, int zF, int yF, int xF, T1 init_value=0) const
#define y0
#define x0
#define xF
#define XMIPP_MIN(x, y)
Definition: xmipp_macros.h:181
#define FINISHINGY(v)
#define yF
#define STARTINGZ(v)

◆ mod()

template<typename T >
void mod ( const MultidimArray< T > &  x,
MultidimArray< T > &  m,
double  y 
)

Definition at line 3902 of file multidim_array.h.

3903 {
3904  m.resizeNoCopy(x);
3905  double *ptr=NULL;
3906  double *ptrm=MULTIDIM_ARRAY(m);
3907  size_t n;
3909  *(ptrm++) = (*ptr) - std::floor((*ptr)/(y))*(y);
3910 }
__host__ __device__ float2 floor(const float2 v)
void resizeNoCopy(const MultidimArray< T1 > &v)
static double * y
#define MULTIDIM_ARRAY(v)
#define FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(v, n, ptr)
int * n

◆ operator!=()

template<typename T >
bool operator!= ( const MultidimArray< T > &  op1,
const MultidimArray< T > &  op2 
)

MultidimArray inequality.

Definition at line 3839 of file multidim_array.h.

3840 {
3841  return !(op1==op2);
3842 }

◆ operator<<() [1/2]

template<typename T >
std::ostream& operator<< ( std::ostream &  ostrm,
const MultidimArray< T > &  v 
)

Output to output stream. This function is not ported to Python.

Definition at line 3916 of file multidim_array.h.

3917 {
3918  if (v.xdim == 0)
3919  ostrm << "NULL Array\n";
3920  else
3921  ostrm << std::endl;
3922 
3923  double max_val = ABS(DIRECT_A3D_ELEM(v , 0, 0, 0));
3924 
3925  T* ptr;
3926  size_t n;
3928  max_val = XMIPP_MAX(max_val, ABS(*ptr));
3929 
3930  int prec = bestPrecision(max_val, 10);
3931 
3932  if (YSIZE(v)==1 && ZSIZE(v)==1)
3933  {
3934  for (int j = STARTINGX(v); j <= FINISHINGX(v); j++)
3935  ostrm << floatToString((double) A3D_ELEM(v, 0, 0, j), 10, prec)
3936  << std::endl;
3937  }
3938  else
3939  {
3940  for (size_t l = 0; l < NSIZE(v); l++)
3941  {
3942  if (NSIZE(v)>1)
3943  ostrm << "Image No. " << l << std::endl;
3944  for (int k = STARTINGZ(v); k <= FINISHINGZ(v); k++)
3945  {
3946  if (ZSIZE(v)>1)
3947  ostrm << "Slice No. " << k << std::endl;
3948  for (int i = STARTINGY(v); i <= FINISHINGY(v); i++)
3949  {
3950  for (int j = STARTINGX(v); j <= FINISHINGX(v); j++)
3951  {
3952  ostrm << floatToString((double) A3D_ELEM(v, k, i, j), 10, prec) << ' ';
3953  }
3954  ostrm << std::endl;
3955  }
3956  }
3957  }
3958  }
3959 
3960  return ostrm;
3961 }
#define NSIZE(v)
#define YSIZE(v)
#define XMIPP_MAX(x, y)
Definition: xmipp_macros.h:193
#define FINISHINGX(v)
int bestPrecision(float F, int _width)
#define FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(v, n, ptr)
#define FINISHINGZ(v)
String floatToString(float F, int _width, int _prec)
#define STARTINGX(v)
#define i
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
#define STARTINGY(v)
#define A3D_ELEM(V, k, i, j)
#define ABS(x)
Definition: xmipp_macros.h:142
#define ZSIZE(v)
#define DIRECT_A3D_ELEM(v, k, i, j)
#define j
#define FINISHINGY(v)
#define STARTINGZ(v)
int * n

◆ operator<<() [2/2]

template<>
std::ostream& operator<< ( std::ostream &  ostrm,
const MultidimArray< std::complex< double > > &  v 
)

Definition at line 86 of file multidim_array.cpp.

88 {
89  if (v.xdim == 0)
90  ostrm << "NULL MultidimArray\n";
91  else
92  ostrm << std::endl;
93 
94  for (size_t l = 0; l < NSIZE(v); l++)
95  {
96  if (NSIZE(v)>1)
97  ostrm << "Image No. " << l << std::endl;
98  for (int k = STARTINGZ(v); k <= FINISHINGZ(v); k++)
99  {
100  if (ZSIZE(v)>1)
101  ostrm << "Slice No. " << k << std::endl;
102  for (int i = STARTINGY(v); i <= FINISHINGY(v); i++)
103  {
104  for (int j = STARTINGX(v); j <= FINISHINGX(v); j++)
105  ostrm << A3D_ELEM(v, k, i, j) << ' ';
106  ostrm << std::endl;
107  }
108  }
109  }
110 
111  return ostrm;
112 }
#define NSIZE(v)
#define FINISHINGX(v)
#define FINISHINGZ(v)
#define STARTINGX(v)
#define i
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
#define STARTINGY(v)
#define A3D_ELEM(V, k, i, j)
#define ZSIZE(v)
#define j
#define FINISHINGY(v)
#define STARTINGZ(v)

◆ operator==() [1/2]

template<typename T >
bool operator== ( const MultidimArray< T > &  op1,
const MultidimArray< T > &  op2 
)

MultidimArray equality.

Definition at line 3832 of file multidim_array.h.

3833 {
3834  return op1.equal(op2);
3835 }
void equal(T op1, MultidimArray< char > &result) const

◆ operator==() [2/2]

template<>
bool operator== ( const MultidimArray< std::complex< double > > &  op1,
const MultidimArray< std::complex< double > > &  op2 
)

Definition at line 162 of file multidim_array.cpp.

163 {
164  double accuracy = XMIPP_EQUAL_ACCURACY;
165  if (! op1.sameShape(op2) || op1.data==NULL || op2.data == NULL)
166  return false;
168  if ( fabs(DIRECT_MULTIDIM_ELEM(op1,n).real() -
169  DIRECT_MULTIDIM_ELEM(op2,n).real() > accuracy)
170  ||
171  fabs(DIRECT_MULTIDIM_ELEM(op1,n).imag() -
172  DIRECT_MULTIDIM_ELEM(op2,n).imag() > accuracy)
173  )
174  return false;
175  return true;
176 }
#define XMIPP_EQUAL_ACCURACY
Definition: xmipp_macros.h:119
bool sameShape(const MultidimArrayBase &op) const
#define FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(v)
#define DIRECT_MULTIDIM_ELEM(v, n)
int * n

◆ planeFit()

void planeFit ( const MultidimArray< double > &  z,
const MultidimArray< double > &  x,
const MultidimArray< double > &  y,
double &  p0,
double &  p1,
double &  p2 
)

Obtains the plane parameters z=p0+p1*x+p2*y.

Definition at line 367 of file multidim_array.cpp.

369 {
371  REPORT_ERROR(ERR_MULTIDIM_SIZE,"Not all vectors are of the same size");
372  if (MULTIDIM_SIZE(z) < 10)
373  REPORT_ERROR(ERR_MULTIDIM_SIZE, "Not enough elements to compute Least Squares plane fit");
374 
375  double m11=0, m12=0, m13=0, m21=0, m22=0, m23=0, m31=0, m32=0, m33=0;
376  double b1=0, b2=0, b3=0;
377 
379  {
380  double X=DIRECT_MULTIDIM_ELEM(x,n);
381  double Y=DIRECT_MULTIDIM_ELEM(y,n);
382  double Z=DIRECT_MULTIDIM_ELEM(z,n);
383  m11+=X*X;
384  m12+=X*Y;
385  m13+=X;
386 
387  m22+=Y*Y;
388  m23+=Y;
389 
390  b1+=X*Z;
391  b2+=Y*Z;
392  b3+=Z;
393  }
394  m21=m12;
395  m31=m13;
396  m32=m23;
397  m33=MULTIDIM_SIZE(z);
398 
399  Matrix2D<double> A(3, 3);
400  Matrix1D<double> b(3);
401  Matrix1D<double> c(3);
402 
403  A(0,0)=m11;
404  A(0,1)=m12;
405  A(0,2)=m13;
406  A(1,0)=m21;
407  A(1,1)=m22;
408  A(1,2)=m23;
409  A(2,0)=m31;
410  A(2,1)=m32;
411  A(2,2)=m33;
412 
413  b(0)=b1;
414  b(1)=b2;
415  b(2)=b3;
416 
417  c = A.inv() * b;
418  p0 = c(2);
419  p2 = c(1);
420  p1 = c(0);
421 }
#define REPORT_ERROR(nerr, ErrormMsg)
Definition: xmipp_error.h:211
doublereal * c
#define MULTIDIM_SIZE(v)
Incorrect MultidimArray size.
Definition: xmipp_error.h:174
doublereal * b
#define FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(v)
#define DIRECT_MULTIDIM_ELEM(v, n)
int * n

◆ selfCoreArrayByArrayMask()

template<typename T >
void selfCoreArrayByArrayMask ( const MultidimArray< T > &  op1,
const MultidimArray< T > &  op2,
MultidimArray< T > &  result,
char  operation,
const MultidimArray< T > *  mask 
)

Core array by array operation.

It assumes that the result is already resized.

Definition at line 2235 of file multidim_array.h.

2238  {
2239  T* ptrResult=NULL;
2240  T* ptrOp1=NULL;
2241  T* ptrOp2=NULL;
2242  T* ptrMask=NULL;
2243  T zero;
2244  zero = T(0);
2245  size_t n;
2246  switch (operation)
2247  {
2248  case '+':
2249  for (n=0, ptrResult=result.data, ptrOp1=op1.data,ptrOp2=op2.data,
2250  ptrMask=mask->data;
2251  n<op1.nzyxdim; n++, ptrResult++,
2252  ptrOp1++, ptrOp2++, ptrMask++)
2253  {
2254  if ((*ptrMask)==zero)
2255  *ptrResult += (*ptrOp1);
2256  else
2257  *ptrResult += (*ptrOp2);//(*ptrOp2) * (*ptrMask);
2258  }
2259  break;
2260  case '-':
2261  for (n=0, ptrResult=result.data, ptrOp1=op1.data,ptrOp2=op2.data,
2262  ptrMask=mask->data;
2263  n<op1.nzyxdim; n++, ptrResult++,
2264  ptrOp1++, ptrOp2++, ptrMask++)
2265  {
2266  if ((*ptrMask)==zero)
2267  *ptrResult -= (*ptrOp1);
2268  else
2269  *ptrResult -= (*ptrOp2) * (*ptrMask);
2270  }
2271  break;
2272  //NOTE: not clear if this should be the default case for * and / operators
2273  case '*':
2274  for (n=0, ptrResult=result.data, ptrOp1=op1.data,ptrOp2=op2.data,
2275  ptrMask=mask->data;
2276  n<op1.nzyxdim; n++, ptrResult++,
2277  ptrOp1++, ptrOp2++, ptrMask++)
2278  {
2279  if ((*ptrMask)==zero)
2280  *ptrResult *= (*ptrOp1);
2281  else
2282  *ptrResult *= (*ptrOp2) * (*ptrMask);
2283  }
2284  break;
2285  case '/':
2286  for (n=0, ptrResult=result.data, ptrOp1=op1.data,ptrOp2=op2.data,
2287  ptrMask=mask->data;
2288  n<op1.nzyxdim; n++, ptrResult++,
2289  ptrOp1++, ptrOp2++, ptrMask++)
2290  {
2291  if ((*ptrMask)==zero)
2292  *ptrResult /= (*ptrOp1);
2293  else
2294  *ptrResult /= (*ptrOp2) * (*ptrMask);
2295  }
2296  break;
2297  }
2298  }
int * n
ql0001_ & zero(ctemp+1),(cvec+1),(a+1),(b+1),(bl+1),(bu+1),(x+1),(w+1), &iout, ifail, &zero,(w+3), &lwar2,(iw+1), &leniw, &glob_grd.epsmac

◆ sincos()

void sincos ( const MultidimArray< double > &  x,
MultidimArray< double > &  s,
MultidimArray< double > &  c 
)

Get Sin and Cos of vector x.

Definition at line 354 of file multidim_array.cpp.

355 {
356  s.resizeNoCopy(x);
357  c.resizeNoCopy(x);
358  double *ptr=NULL;
359  double *ptrS=MULTIDIM_ARRAY(s);
360  double *ptrC=MULTIDIM_ARRAY(c);
361  size_t n;
363  sincos(*ptr, ptrS++,ptrC++);
364 }
void resizeNoCopy(const MultidimArray< T1 > &v)
#define MULTIDIM_ARRAY(v)
#define FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(v, n, ptr)
void sincos(const MultidimArray< double > &x, MultidimArray< double > &s, MultidimArray< double > &c)
int * n

◆ typeCast() [1/3]

template<typename T1 , typename T2 >
void typeCast ( const MultidimArray< T1 > &  v1,
MultidimArray< T2 > &  v2 
)

Conversion from one type to another.

If we have an integer array and we need a double one, we can use this function. The conversion is done through a type casting of each element If n >= 0, only the nth volumes will be converted, otherwise all NSIZE volumes

Definition at line 3749 of file multidim_array.h.

3750 {
3751  if (NZYXSIZE(v1) == 0)
3752  {
3753  v2.clear();
3754  return;
3755  }
3756 
3757  v2.resizeNoCopy(v1);
3758  T1* ptr1 = MULTIDIM_ARRAY(v1);
3759 
3760  // Unroll the loop
3761  const size_t unroll=4;
3762  size_t nmax=(NZYXSIZE(v1)/unroll)*unroll;
3763  for (size_t n=0; n<nmax; n+=unroll, ptr1+=unroll)
3764  {
3765  DIRECT_MULTIDIM_ELEM(v2,n) = static_cast< T2 >(*ptr1);
3766  DIRECT_MULTIDIM_ELEM(v2,n+1) = static_cast< T2 >(*(ptr1+1));
3767  DIRECT_MULTIDIM_ELEM(v2,n+2) = static_cast< T2 >(*(ptr1+2));
3768  DIRECT_MULTIDIM_ELEM(v2,n+3) = static_cast< T2 >(*(ptr1+3));
3769  }
3770  // Do the remaining elements
3771  for (size_t n=nmax; n<NZYXSIZE(v1); ++n, ++ptr1)
3772  DIRECT_MULTIDIM_ELEM(v2,n) = static_cast< T2 >(*ptr1);
3773 }
int * nmax
void resizeNoCopy(const MultidimArray< T1 > &v)
#define MULTIDIM_ARRAY(v)
#define DIRECT_MULTIDIM_ELEM(v, n)
#define NZYXSIZE(v)
int * n

◆ typeCast() [2/3]

template<typename T1 >
void typeCast ( const MultidimArray< T1 > &  v1,
MultidimArray< T1 > &  v2 
)

Conversion from one type to another. In some cases, the two types are the same. So a faster way is simply by assignment.

Definition at line 3814 of file multidim_array.h.

3815 {
3816  v2=v1;
3817 }
double v1

◆ typeCast() [3/3]

template<typename T >
void typeCast ( const MultidimArray< T > &  v1,
Matrix1D< T > &  v2 
)

Assignment.

Definition at line 3822 of file multidim_array.h.

3823 {
3824  v2.resizeNoCopy(XSIZE(v1));
3825  memcpy(&VEC_ELEM(v2,0),&DIRECT_A1D_ELEM(v1,0),MULTIDIM_SIZE(v1)*sizeof(T));
3826  v2.row=false;
3827 }
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
#define MULTIDIM_SIZE(v)
#define DIRECT_A1D_ELEM(v, i)
#define XSIZE(v)
bool row
<0=column vector (default), 1=row vector
Definition: matrix1d.h:267
void resizeNoCopy(int Xdim)
Definition: matrix1d.h:458

◆ typeCastComplex()

template<typename T1 >
void typeCastComplex ( const MultidimArray< T1 > &  v1,
MultidimArray< std::complex< double > > &  v2 
)

Definition at line 3776 of file multidim_array.h.

3777 {
3778  if (NZYXSIZE(v1) == 0)
3779  {
3780  v2.clear();
3781  return;
3782  }
3783 
3784  v2.resizeNoCopy(v1);
3785  T1* ptr1 = MULTIDIM_ARRAY(v1);
3786  double * ptr2 = (double*) MULTIDIM_ARRAY(v2);
3787 
3788  // Unroll the loop
3789  const size_t unroll=4;
3790  size_t nmax=(NZYXSIZE(v1)/unroll)*unroll;
3791  for (size_t n=0; n<nmax; n+=unroll, ptr1+=unroll)
3792  {
3793  *(ptr2++) = static_cast<double>(*ptr1);
3794  *(ptr2++) = 0;
3795  *(ptr2++) = static_cast< double >(*(ptr1+1));
3796  *(ptr2++) = 0;
3797  *(ptr2++) = static_cast< double >(*(ptr1+2));
3798  *(ptr2++) = 0;
3799  *(ptr2++) = static_cast< double >(*(ptr1+3));
3800  *(ptr2++) = 0;
3801  }
3802  // Do the remaining elements
3803  for (size_t n=nmax; n<NZYXSIZE(v1); ++n, ++ptr1)
3804  {
3805  *(ptr2++) = static_cast< double >(*ptr1);
3806  *(ptr2++) = 0;
3807  }
3808 }
int * nmax
void resizeNoCopy(const MultidimArray< T1 > &v)
#define MULTIDIM_ARRAY(v)
#define NZYXSIZE(v)
int * n

◆ window2D()

template<typename T >
void window2D ( const MultidimArray< T > &  Ibig,
MultidimArray< T > &  Ismall,
size_t  y0,
size_t  x0,
size_t  yF,
size_t  xF 
)

Extract piece from image. No check on boundaries are performed.

Definition at line 66 of file multidim_array.cpp.

68 {
69  Ismall.resizeNoCopy(yF - y0 + 1, xF - x0 + 1);
70  STARTINGY(Ismall) = y0;
71  STARTINGX(Ismall) = x0;
72  /*
73  FOR_ALL_ELEMENTS_IN_ARRAY2D(Ismall)
74  A2D_ELEM(Ismall, i, j) = A2D_ELEM(Ibig, i, j);
75  */
76  size_t sizeToCopy=XSIZE(Ismall)*sizeof(T);
77  for (int y=y0; y<=yF; y++)
78  memcpy( &A2D_ELEM(Ismall,y,STARTINGX(Ismall)), &A2D_ELEM(Ibig,y,STARTINGX(Ismall)), sizeToCopy);
79 }
#define A2D_ELEM(v, i, j)
void resizeNoCopy(const MultidimArray< T1 > &v)
static double * y
#define STARTINGX(v)
#define STARTINGY(v)
#define y0
#define x0
#define XSIZE(v)
#define xF
#define yF