Add files/code from https://github.com/manuelruder/artistic-videos

1 files changed, 656 insertions, 0 deletions

diff --git a/consistencyChecker/CTensor4D.h b/consistencyChecker/CTensor4D.h
new file mode 100644
index 0000000..6feeb5d
--- /dev/null
+++ b/consistencyChecker/CTensor4D.h
@@ -0,0 +1,656 @@
+// CTensor4D

+// A four-dimensional array

+//

+// Author: Thomas Brox

+// Last change: 05.11.2001

+//-------------------------------------------------------------------------

+// Note:

+// There is a difference between the GNU Compiler's STL and the standard

+// concerning the definition and usage of string streams as well as substrings.

+// Thus if using a GNU Compiler you should write #define GNU_COMPILER at the

+// beginning of your program.

+//

+// Another Note:

+// Linker problems occured in connection with <vector> from the STL.

+// In this case you should include this file in a namespace.

+// Example:

+// namespace NTensor4D {

+//   #include <CTensor4D.h>

+// }

+// After including other packages you can then write:

+// using namespace NTensor4D;

+

+#ifndef CTENSOR4D_H

+#define CTENSOR4D_H

+

+#include <iostream>

+#include <fstream>

+#include <string>

+#ifdef GNU_COMPILER

+  #include <strstream>

+#else

+  #include <sstream>

+#endif

+#include "CTensor.h"

+

+template <class T>

+class CTensor4D {

+public:

+  // constructor

+  inline CTensor4D();

+  inline CTensor4D(const int aXSize, const int aYSize, const int aZSize, const int aASize);

+  // copy constructor

+  CTensor4D(const CTensor4D<T>& aCopyFrom);

+  // constructor with implicit filling

+  CTensor4D(const int aXSize, const int aYSize, const int aZSize, const int aASize, const T aFillValue);

+  // destructor

+  virtual ~CTensor4D();

+

+  // Changes the size of the tensor, data will be lost

+  void setSize(int aXSize, int aYSize, int aZSize, int aASize);

+  // Downsamples the tensor

+  void downsample(int aNewXSize, int aNewYSize);

+  void downsample(int aNewXSize, int aNewYSize, int aNewZSize);

+  // Upsamples the tensor

+  void upsample(int aNewXSize, int aNewYSize);

+  void upsampleBilinear(int aNewXSize, int aNewYSize);

+  void upsampleTrilinear(int aNewXSize, int aNewYSize, int aNewZSize);

+  // Fills the tensor with the value aValue (see also operator =)

+  void fill(const T aValue);

+  // Copies a box from the tensor into aResult, the size of aResult will be adjusted

+  void cut(CTensor4D<T>& aResult, int x1, int y1, int z1, int a1, int x2, int y2, int z2, int a2);

+  // Reads data from a list of PPM or PGM files given in a text file

+  void readFromFile(char* aFilename);

+  // Writes a set of colour images to a large PPM image

+  void writeToPPM(const char* aFilename, int aCols = 0, int aRows = 0);

+

+  // Gives full access to tensor's values

+  inline T& operator()(const int ax, const int ay, const int az, const int aa) const;

+  // Read access with bilinear interpolation

+  CVector<T> operator()(const float ax, const float ay, const int aa) const;

+  // Fills the tensor with the value aValue (equivalent to fill())

+  inline CTensor4D<T>& operator=(const T aValue);

+  // Copies the tensor aCopyFrom to this tensor (size of tensor might change)

+  CTensor4D<T>& operator=(const CTensor4D<T>& aCopyFrom);

+  // Multiplication with a scalar

+  CTensor4D<T>& operator*=(const T aValue);

+  // Component-wise addition

+  CTensor4D<T>& operator+=(const CTensor4D<T>& aTensor);

+

+  // Gives access to the tensor's size

+  inline int xSize() const;

+  inline int ySize() const;

+  inline int zSize() const;

+  inline int aSize() const;

+  inline int size() const;

+  // Returns the aath layer of the 4D-tensor as 3D-tensor

+  CTensor<T> getTensor3D(const int aa) const;

+  // Removes one dimension and returns the resulting 3D-tensor

+  void getTensor3D(CTensor<T>& aTensor, int aIndex, int aDim = 3) const;

+  // Copies the components of a 3D-tensor in the aDimth layer of the 4D-tensor

+  void putTensor3D(CTensor<T>& aTensor, int aIndex, int aDim = 3);

+    // Removes two dimensions and returns the resulting matrix

+  void getMatrix(CMatrix<T>& aMatrix, int aZIndex, int aAIndex) const;

+  // Copies the components of a 3D-tensor in the aDimth layer of the 4D-tensor

+  void putMatrix(CMatrix<T>& aMatrix, int aZIndex, int aAIndex);

+  // Gives access to the internal data representation (use sparingly)

+  inline T* data() const;

+protected:

+  int mXSize,mYSize,mZSize,mASize;

+  T *mData;

+};

+

+// Provides basic output functionality (only appropriate for very small tensors)

+template <class T> std::ostream& operator<<(std::ostream& aStream, const CTensor4D<T>& aTensor);

+

+// Exceptions thrown by CTensor-------------------------------------------------

+

+// Thrown when one tries to access an element of a tensor which is out of

+// the tensor's bounds

+struct ETensor4DRangeOverflow {

+  ETensor4DRangeOverflow(const int ax, const int ay, const int az, const int aa) {

+    using namespace std;

+    cerr << "Exception ETensor4DRangeOverflow: x = " << ax << ", y = " << ay << ", z = " << az << ", a = " << aa << endl;

+  }

+};

+

+// Thrown from getTensor3D if the parameter's size does not match with the size

+// of this tensor

+struct ETensor4DIncompatibleSize {

+  ETensor4DIncompatibleSize(int ax, int ay, int az, int ax2, int ay2, int az2) {

+    using namespace std;

+    cerr << "Exception ETensor4DIncompatibleSize: x = " << ax << ":" << ax2;

+    cerr << ", y = " << ay << ":" << ay2;

+    cerr << ", z = " << az << ":" << az2 << endl;

+  }

+};

+

+// Thrown from readFromFile if the file format is unknown

+struct ETensor4DInvalidFileFormat {

+  ETensor4DInvalidFileFormat() {

+    using namespace std;

+    cerr << "Exception ETensor4DInvalidFileFormat" << endl;

+  }

+};

+

+// I M P L E M E N T A T I O N --------------------------------------------

+//

+// You might wonder why there is implementation code in a header file.

+// The reason is that not all C++ compilers yet manage separate compilation

+// of templates. Inline functions cannot be compiled separately anyway.

+// So in this case the whole implementation code is added to the header

+// file.

+// Users of CTensor4D should ignore everything that's beyond this line :)

+// ------------------------------------------------------------------------

+

+// P U B L I C ------------------------------------------------------------

+

+// constructor

+template <class T>

+inline CTensor4D<T>::CTensor4D() {

+  mData = 0; mXSize = 0; mYSize = 0; mZSize = 0; mASize = 0;

+}

+

+// constructor

+template <class T>

+inline CTensor4D<T>::CTensor4D(const int aXSize, const int aYSize, const int aZSize, const int aASize)

+  : mXSize(aXSize), mYSize(aYSize), mZSize(aZSize), mASize(aASize) {

+  mData = new T[aXSize*aYSize*aZSize*aASize];

+}

+

+// copy constructor

+template <class T>

+CTensor4D<T>::CTensor4D(const CTensor4D<T>& aCopyFrom)

+  : mXSize(aCopyFrom.mXSize), mYSize(aCopyFrom.mYSize), mZSize(aCopyFrom.mZSize), mASize(aCopyFrom.mASize) {

+  int wholeSize = mXSize*mYSize*mZSize*mASize;

+  mData = new T[wholeSize];

+  for (register int i = 0; i < wholeSize; i++)

+    mData[i] = aCopyFrom.mData[i];

+}

+

+// constructor with implicit filling

+template <class T>

+CTensor4D<T>::CTensor4D(const int aXSize, const int aYSize, const int aZSize, const int aASize, const T aFillValue)

+  : mXSize(aXSize), mYSize(aYSize), mZSize(aZSize), mASize(aASize) {

+  mData = new T[aXSize*aYSize*aZSize*aASize];

+  fill(aFillValue);

+}

+

+// destructor

+template <class T>

+CTensor4D<T>::~CTensor4D() {

+  delete[] mData;

+}

+

+// setSize

+template <class T>

+void CTensor4D<T>::setSize(int aXSize, int aYSize, int aZSize, int aASize) {

+  if (mData != 0) delete[] mData;

+  mData = new T[aXSize*aYSize*aZSize*aASize];

+  mXSize = aXSize;

+  mYSize = aYSize;

+  mZSize = aZSize;

+  mASize = aASize;

+}

+

+//downsample

+template <class T>

+void CTensor4D<T>::downsample(int aNewXSize, int aNewYSize) {

+  T* mData2 = new T[aNewXSize*aNewYSize*mZSize*mASize];

+  int aSize = aNewXSize*aNewYSize;

+  for (int a = 0; a < mASize; a++)

+    for (int z = 0; z < mZSize; z++) {

+      CMatrix<T> aTemp(mXSize,mYSize);

+      getMatrix(aTemp,z,a);

+      aTemp.downsample(aNewXSize,aNewYSize);

+      for (int i = 0; i < aSize; i++)

+        mData2[i+(a*mZSize+z)*aSize] = aTemp.data()[i];

+    }

+  delete[] mData;

+  mData = mData2;

+  mXSize = aNewXSize;

+  mYSize = aNewYSize;

+}

+

+template <class T>

+void CTensor4D<T>::downsample(int aNewXSize, int aNewYSize, int aNewZSize) {

+  T* mData2 = new T[aNewXSize*aNewYSize*aNewZSize*mASize];

+  int aSize = aNewXSize*aNewYSize*aNewZSize;

+  for (int a = 0; a < mASize; a++) {

+    CTensor<T> aTemp(mXSize,mYSize,mZSize);

+    getTensor3D(aTemp,a);

+    aTemp.downsample(aNewXSize,aNewYSize,aNewZSize);

+    for (int i = 0; i < aSize; i++)

+      mData2[i+a*aSize] = aTemp.data()[i];

+  }

+  delete[] mData;

+  mData = mData2;

+  mXSize = aNewXSize;

+  mYSize = aNewYSize;

+  mZSize = aNewZSize;

+}

+

+// upsample

+template <class T>

+void CTensor4D<T>::upsample(int aNewXSize, int aNewYSize) {

+  T* mData2 = new T[aNewXSize*aNewYSize*mZSize*mASize];

+  int aSize = aNewXSize*aNewYSize;

+  for (int a = 0; a < mASize; a++)

+    for (int z = 0; z < mZSize; z++) {

+      CMatrix<T> aTemp(mXSize,mYSize);

+      getMatrix(aTemp,z,a);

+      aTemp.upsample(aNewXSize,aNewYSize);

+      for (int i = 0; i < aSize; i++)

+        mData2[i+(a*mZSize+z)*aSize] = aTemp.data()[i];

+    }

+  delete[] mData;

+  mData = mData2;

+  mXSize = aNewXSize;

+  mYSize = aNewYSize;

+}

+

+// upsampleBilinear

+template <class T>

+void CTensor4D<T>::upsampleBilinear(int aNewXSize, int aNewYSize) {

+  T* mData2 = new T[aNewXSize*aNewYSize*mZSize*mASize];

+  int aSize = aNewXSize*aNewYSize;

+  for (int a = 0; a < mASize; a++)

+    for (int z = 0; z < mZSize; z++) {

+      CMatrix<T> aTemp(mXSize,mYSize);

+      getMatrix(aTemp,z,a);

+      aTemp.upsampleBilinear(aNewXSize,aNewYSize);

+      for (int i = 0; i < aSize; i++)

+        mData2[i+(a*mZSize+z)*aSize] = aTemp.data()[i];

+    }

+  delete[] mData;

+  mData = mData2;

+  mXSize = aNewXSize;

+  mYSize = aNewYSize;

+}

+

+// upsampleTrilinear

+template <class T>

+void CTensor4D<T>::upsampleTrilinear(int aNewXSize, int aNewYSize, int aNewZSize) {

+  T* mData2 = new T[aNewXSize*aNewYSize*aNewZSize*mASize];

+  int aSize = aNewXSize*aNewYSize*aNewZSize;

+  for (int a = 0; a < mASize; a++) {

+    CTensor<T> aTemp(mXSize,mYSize,mZSize);

+    getTensor3D(aTemp,a);

+    aTemp.upsampleTrilinear(aNewXSize,aNewYSize,aNewZSize);

+    for (int i = 0; i < aSize; i++)

+      mData2[i+a*aSize] = aTemp.data()[i];

+  }

+  delete[] mData;

+  mData = mData2;

+  mXSize = aNewXSize;

+  mYSize = aNewYSize;

+  mZSize = aNewZSize;

+}

+

+// fill

+template <class T>

+void CTensor4D<T>::fill(const T aValue) {

+  int wholeSize = mXSize*mYSize*mZSize*mASize;

+  for (register int i = 0; i < wholeSize; i++)

+    mData[i] = aValue;

+}

+

+// cut

+template <class T>

+void CTensor4D<T>::cut(CTensor4D<T>& aResult, int x1, int y1, int z1, int a1, int x2, int y2, int z2, int a2) {

+  aResult.mXSize = x2-x1+1;

+  aResult.mYSize = y2-y1+1;

+  aResult.mZSize = z2-z1+1;

+  aResult.mASize = a2-a1+1;

+  delete[] aResult.mData;

+  aResult.mData = new T[aResult.mXSize*aResult.mYSize*aResult.mZSize*aResult.mASize];

+  for (int a = a1; a <= a2; a++)

+    for (int z = z1; z <= z2; z++)

+      for (int y = y1; y <= y2; y++)

+        for (int x = x1; x <= x2; x++)

+          aResult(x-x1,y-y1,z-z1,a-a1) = operator()(x,y,z,a);

+}

+

+// readFromFile

+template <class T>

+void CTensor4D<T>::readFromFile(char* aFilename) {

+  if (mData != 0) delete[] mData;

+  std::string s;

+  std::string aPath = aFilename;

+  aPath.erase(aPath.find_last_of('\\')+1,100);

+  mASize = 0;

+  {

+    std::ifstream aStream(aFilename);

+    while (!aStream.eof()) {

+      aStream >> s;

+      if (s != "") {

+        mASize++;

+        if (mASize == 1) {

+          s.erase(0,s.find_last_of('.'));

+          if (s == ".ppm" || s == ".PPM") mZSize = 3;

+          else if (s == ".pgm" || s == ".PGM") mZSize = 1;

+          else throw ETensor4DInvalidFileFormat();

+        }

+      }

+    }

+  }

+  std::ifstream aStream(aFilename);

+  aStream >> s;

+  s = aPath+s;

+  // PGM

+  if (mZSize == 1) {

+    CMatrix<float> aTemp;

+    aTemp.readFromPGM(s.c_str());

+    mXSize = aTemp.xSize();

+    mYSize = aTemp.ySize();

+    int aSize = mXSize*mYSize;

+    mData = new T[aSize*mASize];

+    for (int i = 0; i < aSize; i++)

+      mData[i] = aTemp.data()[i];

+    for (int a = 1; a < mASize; a++) {

+      aStream >> s;

+      s = aPath+s;

+      aTemp.readFromPGM(s.c_str());

+      for (int i = 0; i < aSize; i++)

+        mData[i+a*aSize] = aTemp.data()[i];

+    }

+  }

+  // PPM

+  else {

+    CTensor<float> aTemp;

+    aTemp.readFromPPM(s.c_str());

+    mXSize = aTemp.xSize();

+    mYSize = aTemp.ySize();

+    int aSize = 3*mXSize*mYSize;

+    mData = new T[aSize*mASize];

+    for (int i = 0; i < aSize; i++)

+      mData[i] = aTemp.data()[i];

+    for (int a = 1; a < mASize; a++) {

+      aStream >> s;

+      s = aPath+s;

+      aTemp.readFromPPM(s.c_str());

+      for (int i = 0; i < aSize; i++)

+        mData[i+a*aSize] = aTemp.data()[i];

+    }

+  }

+}

+

+// writeToPPM

+template <class T>

+void CTensor4D<T>::writeToPPM(const char* aFilename, int aCols, int aRows) {

+  int rows = (int)floor(sqrt(mASize));

+  if (aRows != 0) rows = aRows;

+  int cols = (int)ceil(mASize*1.0/rows);

+  if (aCols != 0) cols = aCols;

+  FILE* outimage = fopen(aFilename, "wb");

+  fprintf(outimage, "P6 \n");

+  fprintf(outimage, "%ld %ld \n255\n", cols*mXSize,rows*mYSize);

+  for (int r = 0; r < rows; r++)

+    for (int y = 0; y < mYSize; y++)

+      for (int c = 0; c < cols; c++)

+        for (int x = 0; x < mXSize; x++) {

+          unsigned char aHelp;

+          if (r*cols+c >= mASize) aHelp = 0;

+          else aHelp = (unsigned char)operator()(x,y,0,r*cols+c);

+          fwrite (&aHelp, sizeof(unsigned char), 1, outimage);

+          if (r*cols+c >= mASize) aHelp = 0;

+          else aHelp = (unsigned char)operator()(x,y,1,r*cols+c);

+          fwrite (&aHelp, sizeof(unsigned char), 1, outimage);

+          if (r*cols+c >= mASize) aHelp = 0;

+          else aHelp = (unsigned char)operator()(x,y,2,r*cols+c);

+          fwrite (&aHelp, sizeof(unsigned char), 1, outimage);

+        }

+  fclose(outimage);

+}

+

+// operator ()

+template <class T>

+inline T& CTensor4D<T>::operator()(const int ax, const int ay, const int az, const int aa) const {

+  #ifdef DEBUG

+    if (ax >= mXSize || ay >= mYSize || az >= mZSize || aa >= mASize || ax < 0 || ay < 0 || az < 0 || aa < 0)

+      throw ETensorRangeOverflow(ax,ay,az,aa);

+  #endif

+  return mData[mXSize*(mYSize*(mZSize*aa+az)+ay)+ax];

+}

+

+template <class T>

+CVector<T> CTensor4D<T>::operator()(const float ax, const float ay, const int aa) const {

+  CVector<T> aResult(mZSize);

+  int x1 = (int)ax;

+  int y1 = (int)ay;

+  int x2 = x1+1;

+  int y2 = y1+1;

+  #ifdef _DEBUG

+  if (x2 >= mXSize || y2 >= mYSize || x1 < 0 || y1 < 0) throw ETensorRangeOverflow(ax,ay,0);

+  #endif

+  float alphaX = ax-x1; float alphaXTrans = 1.0-alphaX;

+  float alphaY = ay-y1; float alphaYTrans = 1.0-alphaY;

+  for (int k = 0; k < mZSize; k++) {

+    float a = alphaXTrans*operator()(x1,y1,k,aa)+alphaX*operator()(x2,y1,k,aa);

+    float b = alphaXTrans*operator()(x1,y2,k,aa)+alphaX*operator()(x2,y2,k,aa);

+    aResult(k) = alphaYTrans*a+alphaY*b;

+  }

+  return aResult;

+}

+

+// operator =

+template <class T>

+inline CTensor4D<T>& CTensor4D<T>::operator=(const T aValue) {

+  fill(aValue);

+  return *this;

+}

+

+template <class T>

+CTensor4D<T>& CTensor4D<T>::operator=(const CTensor4D<T>& aCopyFrom) {

+  if (this != &aCopyFrom) {

+    if (mData != 0) delete[] mData;

+    mXSize = aCopyFrom.mXSize;

+    mYSize = aCopyFrom.mYSize;

+    mZSize = aCopyFrom.mZSize;

+    mASize = aCopyFrom.mASize;

+    int wholeSize = mXSize*mYSize*mZSize*mASize;

+    mData = new T[wholeSize];

+    for (register int i = 0; i < wholeSize; i++)

+      mData[i] = aCopyFrom.mData[i];

+  }

+  return *this;

+}

+

+// operator *=

+template <class T>

+CTensor4D<T>& CTensor4D<T>::operator*=(const T aValue) {

+  int wholeSize = mXSize*mYSize*mZSize*mASize;

+  for (int i = 0; i < wholeSize; i++)

+    mData[i] *= aValue;

+  return *this;

+}

+

+// operator +=

+template <class T>

+CTensor4D<T>& CTensor4D<T>::operator+=(const CTensor4D<T>& aTensor) {

+  #ifdef _DEBUG

+  if (mXSize != aTensor.mXSize || mYSize != aTensor.mYSize || mZSize != aTensor.mZSize || mASize != aTensor.mASize)

+    throw ETensorIncompatibleSize(mXSize,mYSize,mZSize);

+  #endif

+  int wholeSize = size();

+  for (int i = 0; i < wholeSize; i++)

+    mData[i] += aTensor.mData[i];

+  return *this;

+}

+

+// xSize

+template <class T>

+inline int CTensor4D<T>::xSize() const {

+

+  return mXSize;

+}

+

+// ySize

+template <class T>

+inline int CTensor4D<T>::ySize() const {

+  return mYSize;

+}

+

+// zSize

+template <class T>

+inline int CTensor4D<T>::zSize() const {

+  return mZSize;

+}

+

+// aSize

+template <class T>

+inline int CTensor4D<T>::aSize() const {

+  return mASize;

+}

+

+// size

+template <class T>

+inline int CTensor4D<T>::size() const {

+  return mXSize*mYSize*mZSize*mASize;

+}

+

+// getTensor3D

+template <class T>

+CTensor<T> CTensor4D<T>::getTensor3D(const int aa) const {

+  CTensor<T> aTemp(mXSize,mYSize,mZSize);

+  int aTensorSize = mXSize*mYSize*mZSize;

+  int aOffset = aa*aTensorSize;

+  for (int i = 0; i < aTensorSize; i++)

+    aTemp.data()[i] = mData[i+aOffset];

+  return aTemp;

+}

+

+// getTensor3D

+template <class T>

+void CTensor4D<T>::getTensor3D(CTensor<T>& aTensor, int aIndex, int aDim) const {

+  int aSize;

+  int aOffset;

+  switch (aDim) {

+  case 3:

+    if (aTensor.xSize() != mXSize || aTensor.ySize() != mYSize || aTensor.zSize() != mZSize)

+      throw ETensor4DIncompatibleSize(aTensor.xSize(),aTensor.ySize(),aTensor.zSize(),mXSize,mYSize,mZSize);

+    aSize = mXSize*mYSize*mZSize;

+    aOffset = aIndex*aSize;

+    for (int i = 0; i < aSize; i++)

+      aTensor.data()[i] = mData[i+aOffset];

+    break;

+  case 2:

+    if (aTensor.xSize() != mXSize || aTensor.ySize() != mYSize || aTensor.zSize() != mASize)

+      throw ETensor4DIncompatibleSize(aTensor.xSize(),aTensor.ySize(),aTensor.zSize(),mXSize,mYSize,mASize);

+    aSize = mXSize*mYSize;

+    aOffset = aIndex*aSize;

+    for (int a = 0; a < mASize; a++) 

+      for (int i = 0; i < aSize; i++)

+        aTensor.data()[i+a*aSize] = mData[i+aOffset+a*aSize*mZSize];

+    break;

+  case 1:

+    if (aTensor.xSize() != mXSize || aTensor.ySize() != mZSize || aTensor.zSize() != mASize)

+      throw ETensor4DIncompatibleSize(aTensor.xSize(),aTensor.ySize(),aTensor.zSize(),mXSize,mZSize,mASize);

+    for (int a = 0; a < mASize; a++)

+      for (int z = 0; z < mZSize; z++)

+        for (int x = 0; x < mXSize; x++)

+          aTensor(x,z,a) = operator()(x,aIndex,z,a);

+    break;

+  case 0:

+    if (aTensor.xSize() != mYSize || aTensor.ySize() != mZSize || aTensor.zSize() != mASize)

+      throw ETensor4DIncompatibleSize(aTensor.xSize(),aTensor.ySize(),aTensor.zSize(),mYSize,mZSize,mASize);

+    for (int a = 0; a < mASize; a++)

+      for (int z = 0; z < mZSize; z++)

+        for (int y = 0; y < mYSize; y++)

+          aTensor(y,z,a) = operator()(aIndex,y,z,a);

+    break;

+  default: getTensor3D(aTensor,aIndex);

+  }

+}

+

+// putTensor3D

+template <class T>

+void CTensor4D<T>::putTensor3D(CTensor<T>& aTensor, int aIndex, int aDim) {

+  int aSize;

+  int aOffset;

+  switch (aDim) {

+  case 3:

+    if (aTensor.xSize() != mXSize || aTensor.ySize() != mYSize || aTensor.zSize() != mZSize)

+      throw ETensor4DIncompatibleSize(aTensor.xSize(),aTensor.ySize(),aTensor.zSize(),mXSize,mYSize,mZSize);

+    aSize = mXSize*mYSize*mZSize;

+    aOffset = aIndex*aSize;

+    for (int i = 0; i < aSize; i++)

+      mData[i+aOffset] = aTensor.data()[i];

+    break;

+  case 2:

+    if (aTensor.xSize() != mXSize || aTensor.ySize() != mYSize || aTensor.zSize() != mASize)

+      throw ETensor4DIncompatibleSize(aTensor.xSize(),aTensor.ySize(),aTensor.zSize(),mXSize,mYSize,mASize);

+    aSize = mXSize*mYSize;

+    aOffset = aIndex*aSize;

+    for (int a = 0; a < mASize; a++)

+      for (int i = 0; i < aSize; i++)

+        mData[i+aOffset+a*aSize*mZSize] = aTensor.data()[i+a*aSize];

+    break;

+  case 1:

+    if (aTensor.xSize() != mXSize || aTensor.ySize() != mZSize || aTensor.zSize() != mASize)

+      throw ETensor4DIncompatibleSize(aTensor.xSize(),aTensor.ySize(),aTensor.zSize(),mXSize,mZSize,mASize);

+    for (int a = 0; a < mASize; a++)

+      for (int z = 0; z < mZSize; z++)

+        for (int x = 0; x < mXSize; x++)

+          operator()(x,aIndex,z,a) = aTensor(x,z,a);

+    break;

+  case 0:

+    if (aTensor.xSize() != mYSize || aTensor.ySize() != mZSize || aTensor.zSize() != mASize)

+      throw ETensor4DIncompatibleSize(aTensor.xSize(),aTensor.ySize(),aTensor.zSize(),mYSize,mZSize,mASize);

+    for (int a = 0; a < mASize; a++)

+      for (int z = 0; z < mZSize; z++)

+        for (int y = 0; y < mYSize; y++)

+          operator()(aIndex,y,z,a) = aTensor(y,z,a);

+    break;

+  default: putTensor3D(aTensor,aIndex);

+  }

+}

+

+// getMatrix

+template <class T>

+void CTensor4D<T>::getMatrix(CMatrix<T>& aMatrix, int aZIndex, int aAIndex) const {

+  if (aMatrix.xSize() != mXSize || aMatrix.ySize() != mYSize)

+    throw ETensor4DIncompatibleSize(aMatrix.xSize(),aMatrix.ySize(),1,mXSize,mYSize,1);

+  int aSize = mXSize*mYSize;

+  int aOffset = aSize*(aAIndex*mZSize+aZIndex);

+  for (int i = 0; i < aSize; i++)

+    aMatrix.data()[i] = mData[i+aOffset];

+}

+

+// putMatrix

+template <class T>

+void CTensor4D<T>::putMatrix(CMatrix<T>& aMatrix, int aZIndex, int aAIndex) {

+  if (aMatrix.xSize() != mXSize || aMatrix.ySize() != mYSize)

+    throw ETensor4DIncompatibleSize(aMatrix.xSize(),aMatrix.ySize(),1,mXSize,mYSize,1);

+  int aSize = mXSize*mYSize;

+  int aOffset = aSize*(aAIndex*mZSize+aZIndex);

+  for (int i = 0; i < aSize; i++)

+    mData[i+aOffset] = aMatrix.data()[i];

+}

+

+// data()

+template <class T>

+inline T* CTensor4D<T>::data() const {

+  return mData;

+}

+

+// N O N - M E M B E R  F U N C T I O N S --------------------------------------

+

+// operator <<

+template <class T>

+std::ostream& operator<<(std::ostream& aStream, const CTensor4D<T>& aTensor) {

+  for (int a = 0; a < aTensor.aSize(); a++) {

+    for (int z = 0; z < aTensor.zSize(); z++) {

+      for (int y = 0; y < aTensor.ySize(); y++) {

+        for (int x = 0; x < aTensor.xSize(); x++)

+          aStream << aTensor(x,y,z) << ' ';

+        aStream << std::endl;

+      }

+      aStream << std::endl;

+    }

+    aStream << std::endl;

+  }

+  return aStream;

+}

+

+#endif