merged cmake with master

core/basics/binstream.h

@@ -325,11 +325,12 @@ inline obinstream& operator<<(obinstream &s, const LinAl::MatrixC<Tp> &A) {
   return s;
 }
 
-#else
-#ifndef LINAL_WARNING
-#pragma message NICE_WARNING("LinAl addons will not be compiled.")
-#define LINAL_WARNING
-#endif
+//#else
+//#ifndef LINAL_WARNING
+//#pragma message NICE_WARNING("LinAl addons will not be compiled.")
+//#define LINAL_WARNING
+//#endif
+
 #endif
 
  

core/image/MultiChannelImage3DT.cpp

@@ -0,0 +1,34 @@
+#include "core/image/MultiChannelImage3DT.h"
+
+using namespace NICE;
+
+#define SET_FUNCS_MACRO(MYTYPE) \
+  template <> \
+  int MultiChannelImage3DT<MYTYPE>::getPixelInt(int x, int y, int z, int channel) const \
+  { \
+    return (int)get(x,y,z,channel); \
+  } \
+  \
+  template<>\
+  double MultiChannelImage3DT<MYTYPE>::getPixelFloat(int x, int y, int z, int channel) const\
+  {\
+    return (double)get(x,y,z,channel);\
+  }\
+  template<>\
+  void MultiChannelImage3DT<MYTYPE>::setPixelInt(int x, int y, int z, int channel, int pixel)\
+  {\
+    set(x,y,z,(MYTYPE)pixel,channel);\
+  }\
+  \
+  template<>\
+  void MultiChannelImage3DT<MYTYPE>::setPixelFloat(int x, int y, int z, int channel, double pixel)\
+  {\
+    set(x,y,z,(MYTYPE)pixel,channel);\
+  }
+
+SET_FUNCS_MACRO ( int )
+SET_FUNCS_MACRO ( long int )
+SET_FUNCS_MACRO ( double )
+SET_FUNCS_MACRO ( float )
+SET_FUNCS_MACRO ( unsigned int )
+

core/image/MultiChannelImage3DT.h

@@ -0,0 +1,183 @@
+#ifndef _NICE_CORE_MULTICHANNELIMAGE3DT_H
+#define _NICE_CORE_MULTICHANNELIMAGE3DT_H
+
+#include <core/image/MultiChannelImageAccess3D.h>
+#include <core/image/ImageT.h>
+#include <core/image/MultiChannelImageT.h>
+
+#include <vector>
+#include <fstream>
+
+namespace NICE {
+
+/**
+ * @class MultiChannelImage3DT
+ * A 3d image (arbitrary number of cross section images) consisting of an arbitrary number of channels.
+ *
+ * formaly known as Generic Image
+ *
+ * @author Björn Fröhlich and Sven Sickert
+ * @example notyet
+ */
+
+template <class P>
+class MultiChannelImage3DT : public MultiChannelImageAccess3D {
+protected:
+  typedef P Value;
+  typedef unsigned int uint;
+
+  /** image data, use carefully !!! data[channel][pixel_offset] */
+  P **data;
+  
+  /** image width */
+  int xsize;
+  
+  /** image height */
+  int ysize;
+
+  /** image depth */
+  int zsize;
+  
+  /** number of image channels */
+  uint numChannels;
+  
+public:
+
+  virtual inline int width() const;
+  virtual inline int height() const;
+  virtual inline int depth() const;
+
+  virtual inline int channels() const;
+
+  virtual int getPixelInt( int x, int y, int z, int channel ) const;
+
+  virtual double getPixelFloat( int x, int y, int z, int channel ) const;
+
+  virtual void setPixelInt( int x, int y, int z, int channel, int pixel );
+
+  virtual void setPixelFloat( int x, int y, int z, int channel, double pixel );
+
+  /** simple constructor */
+  MultiChannelImage3DT( int xsize, int ysize, int zsize, uint numChannels = 1);
+
+  /** very simple constructor */
+  MultiChannelImage3DT();
+
+  /** copy constructor */
+  MultiChannelImage3DT( const MultiChannelImage3DT<P>& p );
+
+  /** simple destructor */
+  virtual ~MultiChannelImage3DT();
+
+  /** free all memory */
+  void freeData();
+
+  /** reinit */
+  void reInit( int xsize, int ysize, int zsize, int numChannels = 1);
+
+  /** reinit data structure using the same dimensions and
+   number of channels as another image */
+  template<class SrcP>
+  void reInitFrom( const MultiChannelImage3DT<SrcP> & src);
+
+  void addChannel( int newChans = 1 );
+
+  /** add a channel to Multichannel Image */
+  template<class SrcP>
+  void addChannel( NICE::MultiChannelImageT<SrcP> &newImg );
+  
+  template<class SrcP>
+  void addChannel(const NICE::MultiChannelImage3DT<SrcP> &newImg);
+  
+  /** get value */
+  P get( int x, int y, int z, uint channel = 0 ) const;
+  
+  /** get data pointer */
+  P** getDataPointer();
+
+  /** set value */
+  void set( int x, int y, int z, P val, uint channel = 0 );
+
+  /** set value */
+  void set( P val, uint channel = 0 );
+
+  /** set value */
+  void setAll( P val );
+
+  /** calc integral image */
+  void calcIntegral( uint channel = 0 );
+  
+  /** 
+   * @brief calculate the gray level co-occurence matrix of a channel
+   * @param max Co-occurence Matrix
+   * @param dis displacement vector
+   * @param channel channel
+   */
+  void calcGLCM( std::vector<std::vector<double> > & mat, const std::vector<int> dis, uint channel = 0 );
+  
+  /** 
+   * @brief calculate the variance image map of a channel
+   * @param srcchan source channel with raw data
+   * @param tarchan target channel for the variance map
+   */
+  void calcVariance( uint srcchan = 0, uint tarchan = 1 );
+  
+  /**
+   * @brief calculate the integral value in the volume given by upper left front corner and lower right back corner, including out of boundary check
+   * @warning make sure that the given channel is an integral 3d image
+   * @param ulfx upper left front x coordinate
+   * @param ulfy upper left front y coordinate
+   * @param ulfz upper left front z coordinate
+   * @param lrbx lower right back x coordinate
+   * @param lrby lower right back y coordinate
+   * @param lrbz lower right back z coordinate
+   * @param channel channel
+   * @return P mean value of given volume
+   **/
+  P getIntegralValue(int ulfx, int ulfy, int ulfz, int lrbx, int lrby, int lrbz, int channel);
+
+  /** convert to ice image */
+  void convertToGrey( NICE::Image & img, int z, uint channel = 0, bool normalize = true ) const;
+
+	/** convert to ice image template */
+  void convertToGrey( NICE::ImageT<P> & img, int z, uint channel = 0, bool normalize = false ) const;
+	
+  /** convert to ice colorimage */
+  void convertToColor( NICE::ColorImage & img, int z, const int chan1 = 0, const int chan2 = 1, const int chan3 = 2 ) const;
+
+  /** return image for visualization */
+  Image getChannel( int z, uint channel = 0 ) const;
+
+	/** return image for visualization */
+  ImageT<P> getChannelT( int z, uint channel = 0 ) const;
+
+  /** return rgb image for visualization */
+  ColorImage getColor(int z) const;
+
+  /** calculate image statistics */
+  void statistics( P & min, P & max, uint channel = 0 ) const;
+
+  /** correct inhomogeneous illuminations (shading) between the image slices **/
+  void correctShading( uint channel = 0 ) const;
+  
+  /** dump all data to RAW format: xsize, ysize, numChannels, <data> */
+  void store( std::string filename ) const;
+
+  /** read all data from RAW format: xsize, ysize, numChannels, <data> */
+  void restore( std::string filename );
+  
+  /** copy alls data to new object */
+  MultiChannelImage3DT<P>& operator=( const MultiChannelImage3DT<P>& orig );
+  
+  /** element operator */
+  P & operator() (int x, int y, int z, uint channel = 0);
+  
+  /** element operator */
+  MultiChannelImageT<P> operator[] (uint c);
+};
+
+} // namespace
+
+#include <core/image/MultiChannelImage3DT.tcc>
+
+#endif

core/image/MultiChannelImage3DT.tcc

@@ -0,0 +1,920 @@
+#include <iostream>
+#include <assert.h>
+#include <stdio.h>
+#include <vector>
+
+namespace NICE {
+template<class P>
+MultiChannelImage3DT<P>::MultiChannelImage3DT( int _xsize, int _ysize, int _zsize, uint _numChannels)
+{
+  data = NULL;
+  numChannels = 0;
+  xsize = 0;
+  ysize = 0;
+  zsize = 0;
+  reInit( _xsize, _ysize, _zsize, _numChannels);
+}
+
+template<class P>
+MultiChannelImage3DT<P>::MultiChannelImage3DT()
+{
+  xsize = 0;
+  ysize = 0;
+  zsize = 0;
+  numChannels = 0;
+  data = NULL;
+}
+
+template<class P>
+P & MultiChannelImage3DT<P>::operator() (int x, int y, int z, uint channel)
+{
+  assert( channel < numChannels );
+  assert(( x < xsize ) && ( x >= 0 ) );
+  assert(( y < ysize ) && ( y >= 0 ) );
+  assert(( z < zsize ) && ( z >= 0 ) );
+  assert( data[channel] != NULL );
+  return data[channel][x + y*xsize + z*xsize*ysize];
+}
+
+template<class P>
+MultiChannelImageT<P> MultiChannelImage3DT<P>::operator[] (uint c)
+{
+  // This was our first idea ... but it creates a real copy
+  // ImageT<P> tmp ( data[c], xsize, ysize, xsize*sizeof(P), GrayColorImageCommonImplementation::noAlignment );
+  // This is the correct version. The funny thing about this is that shallowCopy
+  // is not an enum parameter, but an instance of ShallowCopyMode, which is a class.
+  // This fancy trick was done in older to prevent automatic conversion between enum types
+  // as done implicitly by C++.
+  
+  MultiChannelImageT<P> img;
+  for( int z = 0; z < zsize; z++ )
+  {
+    P * datatmp = data[c];
+    ImageT<P> tmp ( &datatmp[z*(xsize*ysize)], xsize, ysize, xsize*sizeof(P), GrayColorImageCommonImplementation::shallowCopy );
+    img.addChannel(tmp);
+  }
+  return img;
+}
+
+template<class P>
+MultiChannelImage3DT<P>& MultiChannelImage3DT<P>::operator=(const MultiChannelImage3DT<P>& orig) 
+{
+  if(!(xsize == orig.xsize && ysize == orig.ysize && zsize == orig.zsize && numChannels == orig.numChannels))
+  {
+    freeData();
+    xsize = orig.xsize;
+    ysize = orig.ysize;
+    zsize = orig.zsize;
+    numChannels = orig.numChannels;
+    if(orig.data != NULL)
+    {
+      data = new P *[numChannels];
+      for ( int c = 0; c < ( int )numChannels; c++ )
+      {
+        if ( orig.data[c] == NULL )
+        {
+          data[c] = NULL;
+        }
+        else
+        {
+          data[c] = new P [xsize*ysize*zsize];
+        }
+      }
+    }
+    else
+      data = NULL;
+  }
+
+  for ( int c = 0; c < ( int )numChannels; c++ )
+  {
+    if ( orig.data[c] != NULL )
+    {
+      for ( int x = 0; x < xsize*ysize*zsize; x++ )
+      {
+        data[c][x] = orig.data[c][x];
+      }
+    }
+  }
+  
+  return *this;
+}
+
+template<class P>
+MultiChannelImage3DT<P>::MultiChannelImage3DT( const MultiChannelImage3DT<P>& p )
+{
+  data = NULL;
+  xsize = p.xsize;
+  ysize = p.ysize;
+  zsize = p.zsize;
+  numChannels = p.numChannels;
+   
+  if(p.data != NULL)
+    data = new P *[numChannels];
+  else
+    data = NULL;
+
+  for ( int c = 0; c < ( int )numChannels; c++ )
+  {
+    if ( p.data[c] == NULL )
+    {
+      data[c] = NULL;
+    }
+    else
+    {
+      data[c] = new P [xsize*ysize*zsize];
+      for ( int x = 0; x < xsize*ysize*zsize; x++ )
+      {
+        data[c][x] = p.data[c][x];
+      }
+    }
+  }
+}
+
+
+template<class P>
+void MultiChannelImage3DT<P>::addChannel( int newChans )
+{
+  P **tmpData = new P *[numChannels+newChans];
+
+  bool allocMem = false;
+  int i = 0;
+
+  for ( ; i < (int)numChannels; i++ )
+  {
+    tmpData[i] = data[i];
+
+    if ( data[i] != NULL )
+      allocMem = true;
+  }
+
+  if ( allocMem )
+  {
+    for ( ; i < newChans + (int)numChannels; i++ )
+    {
+      tmpData[i] = new P [xsize*ysize*zsize];
+    }
+  }
+
+  numChannels += newChans;
+
+  delete [] data;
+
+  data = new P *[numChannels];
+
+  for ( i = 0; i < (int)numChannels; i++ )
+  {
+    data[i] = tmpData[i];
+  }
+
+  delete [] tmpData;
+}
+
+template<class P>
+template<class SrcP>
+void MultiChannelImage3DT<P>::addChannel(NICE::MultiChannelImageT<SrcP> &newMCImg)
+{
+  int oldchan = numChannels;
+  if(this->xsize > 0)
+  {
+    assert(newMCImg.width() == this->width() && newMCImg.height() == this->height());
+    assert(newMCImg.channels() == this->zsize);
+    addChannel(1);
+  }
+  else
+  {
+    reInit( newMCImg.width(), newMCImg.height(), newMCImg.channels(), 1 );
+  }
+  
+  for(int z = 0; z < this->zsize; z++)
+  {
+    NICE::ImageT<SrcP> newImg = newMCImg[z];
+    for(int y = 0; y < this->ysize; y++)
+    {
+      for(int x = 0; x < this->xsize; x++)
+      {
+        data[oldchan][x + y*xsize + z*xsize*ysize] = (P)newImg(x,y);
+      }
+    }
+  }
+}
+
+template<class P>
+template<class SrcP>
+void MultiChannelImage3DT<P>::addChannel(const NICE::MultiChannelImage3DT<SrcP> &newImg)
+{
+  int oldchan = numChannels;
+  if(numChannels > 0)
+  {
+    assert(newImg.width() == this->width() && newImg.height() == this->height() && newImg.depth() == this->depth());
+    addChannel(newImg.channels());
+  }
+  else
+  {
+      reInit( newImg.width(), newImg.height(), newImg.depth(), newImg.channels() );
+  }
+  
+  int chanNI = 0;
+  
+  for(int c = oldchan; c < (int)numChannels; c++, chanNI++)
+  {
+    int val = 0;
+    for(int z = 0; z < this->zsize; z++)
+    {
+      for(int y = 0; y < this->ysize; y++)
+      {
+        for(int x = 0; x < this->xsize; x++, val++)
+        {
+          data[c][val] = newImg.get(x,y,z,chanNI);
+        }
+      }
+    }
+  }
+}
+
+template<class P>
+MultiChannelImage3DT<P>::~MultiChannelImage3DT()
+{
+  freeData();
+}
+
+template<class P>
+void MultiChannelImage3DT<P>::freeData()
+{
+  if ( data != NULL )
+  {
+    for ( uint i = 0 ; i < numChannels ; i++ )
+      if ( data[i] != NULL )
+        delete [] data[i];
+
+    delete [] data;
+
+    data = NULL;
+  }
+}
+
+template<class P>
+void MultiChannelImage3DT<P>::reInit( int _xsize, int _ysize, int _zsize, int _numChannels )
+{
+  freeData();
+  xsize = _xsize;
+  ysize = _ysize;
+  zsize = _zsize;
+  numChannels = _numChannels;
+  data = new P *[numChannels];
+
+  for ( uint i = 0 ; i < numChannels; i++ )
+    data[i] = new P [xsize*ysize*zsize];
+}
+
+template<class P>
+template<class SrcP>
+void MultiChannelImage3DT<P>::reInitFrom( const MultiChannelImage3DT<SrcP> & src )
+{
+  freeData();
+
+  xsize = src.width();
+  ysize = src.height();
+  zsize = src.depth();
+  numChannels = src.channels();
+
+  data = new P *[numChannels];
+
+  for ( uint i = 0 ; i < numChannels; i++ )
+    data[i] = new P [xsize*ysize*zsize];
+}
+
+template<class P>
+P MultiChannelImage3DT<P>::get( int x, int y, int z, uint channel ) const
+{
+  assert( channel < numChannels );
+  assert(( x < xsize ) && ( x >= 0 ) );
+  assert(( y < ysize ) && ( y >= 0 ) );
+  assert(( z < zsize ) && ( z >= 0 ) );
+  assert( data[channel] != NULL );
+
+  return data[channel][x + y*xsize + z*xsize*ysize];
+}
+
+template<class P>
+P ** MultiChannelImage3DT<P>::getDataPointer()
+{
+  return data;
+}
+
+template<class P>
+void MultiChannelImage3DT<P>::set( int x, int y, int z, P val, uint channel )
+{
+  assert( channel < numChannels );
+  assert(( x < xsize ) && ( x >= 0 ) );
+  assert(( y < ysize ) && ( y >= 0 ) );
+  assert(( z < zsize ) && ( z >= 0 ) );
+  assert( data[channel] != NULL );
+
+  data[channel][x + y*xsize + z*xsize*ysize] = val;
+}
+
+template<class P>
+void MultiChannelImage3DT<P>::set( P val, uint channel )
+{
+  assert( channel < numChannels );
+  assert( data[channel] != NULL );
+
+  for ( int k = 0 ; k < xsize*ysize*zsize ; k++ )
+    data[channel][k] = val;
+}
+
+template<class P>
+void MultiChannelImage3DT<P>::setAll( P val )
+{
+  for ( uint channel = 0 ; channel < numChannels ; channel++ )
+    if ( data[channel] != NULL )
+      set( val, channel );
+}
+
+template<class P>
+void MultiChannelImage3DT<P>::statistics( P & min, P & max, uint channel ) const
+{
+  assert( channel < numChannels );
+
+  for ( long k = 0 ; k < xsize*ysize*zsize ; k++ )
+  {
+    P val = data [channel][k];
+    if (( k == 0 ) || ( val > max ) ) max = val;
+    if (( k == 0 ) || ( val < min ) ) min = val;
+  }
+  
+  assert(finite(max));
+  assert(finite(min));
+}
+
+template<class P>
+void MultiChannelImage3DT<P>::correctShading( uint channel ) const
+{
+  assert( channel < numChannels );
+  
+  std::vector<double> meanVals;
+  for( int z = 0; z < zsize; z++ )
+  {
+    double sumVal = 0;
+    for( int y = 0; y < ysize; y++ )
+    {
+      for( int x = 0; x < xsize; x++ )
+      {
+        sumVal += data [channel][x + y*xsize + z*xsize*ysize];
+      }
+    }
+    sumVal /= (xsize*ysize);
+    meanVals.push_back( sumVal );
+  }
+
+  P newMax = std::numeric_limits<P>::min();
+  short int maxVal = 255;
+  for ( int z = 0; z < zsize; z++ )
+  {
+    for ( int y = 0; y < ysize; y++ )
+    {
+      for ( int x = 0; x < xsize; x++ )
+      {
+        P tmp = data [channel][x + y*xsize + z*xsize*ysize];
+        double newVal = maxVal * ( (double) tmp / meanVals[z] );
+        if ( ( P ) newVal > newMax )
+          newMax = ( P ) newVal;
+        data [channel][x + y*xsize + z*xsize*ysize] = newVal;
+      }
+    }
+  }
+  
+  for ( long k = 0 ; k < xsize*ysize*zsize ; k++ )
+  {
+    data [channel][k] = data [channel][k] / newMax * maxVal;
+  }
+}
+
+template<class P>
+Image MultiChannelImage3DT<P>::getChannel( int z, uint channel ) const
+{
+  assert( channel < numChannels );
+
+  NICE::Image img;
+  convertToGrey( img, z, channel, true );
+
+  return img;
+}
+
+template<class P>
+ImageT<P> MultiChannelImage3DT<P>::getChannelT( int z, uint channel ) const
+{
+  assert( channel < numChannels );
+
+  NICE::ImageT<P> img;
+  convertToGrey( img, z, channel, false );
+  
+  P min, max;
+  statistics ( min, max, channel );
+  fprintf (stderr, "MultiChannelImage3DT<>::showChannel: max %f min %f\n", (double)max, (double)min );
+
+  return img;
+}
+
+/** convert to ice image */
+template<class P>
+void MultiChannelImage3DT<P>::convertToGrey( NICE::Image & img, int z, uint channel,  bool normalize ) const
+{
+  assert( channel < numChannels );
+
+  P min, max;
+
+  if ( normalize ) {
+    statistics( min, max, channel );
+    fprintf( stderr, "MultiChannelImage3DT<>::showChannel: max %f min %f\n", ( double )max, ( double )min );
+  }
+
+  bool skip_assignment = false;
+
+  img.resize( xsize, ysize );
+
+  if ( normalize )
+    if ( max - min < std::numeric_limits<double>::min() )
+    {
+      img.set( max );
+      skip_assignment = true;
+      fprintf( stderr, "MultiChannelImage3DT::showChannel: image is uniform! (%f)\n", ( double )max );
+    }
+
+
+  if ( ! skip_assignment )
+  {
+    long k = 0;
+
+    for ( int y = 0 ; y < ysize; y++ )
+    {
+      for ( int x = 0 ; x < xsize ; x++, k++ )
+      {
+        if ( normalize )
+        {
+          img.setPixel( x, y, ( int )(( data[channel][z*xsize*ysize + k] - min ) * 255 / ( max - min ) ) );
+        }
+        else
+        {
+          img.setPixel( x, y, ( int )( data[channel][z*xsize*ysize + k] ) );
+        }
+      }
+    }
+  }
+}
+
+/** convert to ice image template */
+template<class P>
+void MultiChannelImage3DT<P>::convertToGrey( NICE::ImageT<P> & img, int z, uint channel,  bool normalize ) const
+{
+  assert( channel < numChannels );
+
+  P min, max;
+
+  if ( normalize ) {
+    statistics( min, max, channel );
+    fprintf( stderr, "MultiChannelImage3DT<>::showChannel: max %f min %f\n", ( double )max, ( double )min );
+  }
+
+  bool skip_assignment = false;
+
+  img.resize( xsize, ysize );
+
+  if ( normalize )
+    if ( max - min < std::numeric_limits<double>::min() )
+    {
+      img.set( max );
+      skip_assignment = true;
+      fprintf( stderr, "MultiChannelImage3DT::showChannel: image is uniform! (%f)\n", ( double )max );
+    }
+
+
+  if ( ! skip_assignment )
+  {
+    long k = 0;
+
+    for ( int y = 0 ; y < ysize; y++ )
+    {
+      for ( int x = 0 ; x < xsize ; x++, k++ )
+      {
+        if ( normalize )
+        {
+          img.setPixel( x, y, ( int )(( data[channel][z*xsize*ysize + k] - min ) * 255 / ( max - min ) ) );
+        }
+        else
+        {
+          img.setPixel( x, y, ( int )( data[channel][z*xsize*ysize + k] ) );
+        }
+      }
+    }
+  }
+}
+
+template<class P>
+void MultiChannelImage3DT<P>::convertToColor( NICE::ColorImage & img, int z, const int chan1, const int chan2, const int chan3) const
+{
+  assert( chan1 < numChannels && chan2 < numChannels && chan3 < numChannels);
+
+  img.resize( xsize, ysize );
+
+  long k = 0;
+
+  for ( int y = 0 ; y < ysize; y++ )
+  {
+    for ( int x = 0 ; x < xsize ; x++, k++ )
+    {
+      img.setPixel( x, y, 0, ( int )( data[chan1][z*xsize*ysize + k] ) );
+      img.setPixel( x, y, 1, ( int )( data[chan2][z*xsize*ysize + k] ) );
+      img.setPixel( x, y, 2, ( int )( data[chan3][z*xsize*ysize + k] ) );
+    }
+  }
+}
+
+template<class P>
+ColorImage MultiChannelImage3DT<P>::getColor(int z) const
+{
+  assert( z < zsize );
+	assert( numChannels >= 3 );
+
+  NICE::ColorImage img( xsize, ysize );
+
+  long k = 0;
+
+  for ( int y = 0 ; y < ysize; y++ )
+  {
+    for ( int x = 0 ; x < xsize ; x++, k++ )
+    {
+      img.setPixel( x, y, 0, ( int )( data[0][z*xsize*ysize + k] ) );
+      img.setPixel( x, y, 1, ( int )( data[1][z*xsize*ysize + k] ) );
+      img.setPixel( x, y, 2, ( int )( data[2][z*xsize*ysize + k] ) );
+    }
+  }
+  //showImage(img);
+  //getchar();
+  return img;
+}
+
+template<class P>
+void MultiChannelImage3DT<P>::calcIntegral( uint channel )
+{
+  assert( channel < numChannels );
+  assert( data[channel] != NULL );
+
+  P *integralImage = data[channel];
+
+  /** first column **/
+  int k = xsize;
+  for ( int y = 1 ; y < ysize; y++, k += xsize )
+    integralImage[k] += integralImage[k-xsize];
+
+  /** first row **/
+  k = 1;
+  for ( int x = 1 ; x < xsize; x++, k++ )
+    integralImage[k] += integralImage[k-1];
+
+  /** first stack (depth) **/
+  k = xsize * ysize;
+  for ( int z = 1 ; z < zsize; z++, k += (xsize*ysize) )
+    integralImage[k] += integralImage[k-(xsize*ysize)];
+
+  /** x-y plane **/
+  k = xsize + 1;
+  for ( int y = 1 ; y < ysize ; y++, k++ )
+    for ( int x = 1 ; x < xsize ; x++, k++ )
+    {
+      integralImage[k] += integralImage[k-1];
+      integralImage[k] += integralImage[k - xsize];
+      integralImage[k] -= integralImage[k - xsize - 1];
+    }
+
+  /** y-z plane **/
+  k = xsize*ysize + xsize;
+  for ( int z = 1 ; z < zsize ; z++, k+=xsize )
+    for ( int y = 1 ; y < zsize ; y++, k+=xsize )
+    {
+      integralImage[k] += integralImage[k-(xsize*ysize)];
+      integralImage[k] += integralImage[k - xsize];
+      integralImage[k] -= integralImage[k - xsize - (xsize*ysize)];
+    }
+  
+  /** x-z plane **/
+  k = xsize*ysize + 1;
+  for ( int z = 1 ; z < zsize ; z++, k+=((xsize*ysize)-(xsize-1)) )
+    for ( int x = 1 ; x < xsize ; x++, k++ )
+    {
+      integralImage[k] += integralImage[k-1];
+      integralImage[k] += integralImage[k - (xsize*ysize)];
+      integralImage[k] -= integralImage[k - (xsize*ysize) - 1];
+    }
+    
+  /** all other pixels **/
+  k = xsize*ysize + xsize + 1;
+  for ( int z = 1 ; z < zsize ; z++, k+= xsize )
+  {
+    for ( int y = 1 ; y < ysize ; y++, k++ )
+    {
+      for ( int x = 1 ; x < xsize ; x++, k++ )
+      {
+        integralImage[k] += integralImage[k - (xsize*ysize)];
+        integralImage[k] += integralImage[k - xsize];
+        integralImage[k] += integralImage[k - 1];
+        integralImage[k] += integralImage[k - (xsize*ysize) - xsize - 1];
+        integralImage[k] -= integralImage[k - (xsize*ysize) - xsize];
+        integralImage[k] -= integralImage[k - (xsize*ysize) - 1];
+        integralImage[k] -= integralImage[k - xsize - 1];
+      }
+    }
+  }
+}
+
+template<class P>
+void MultiChannelImage3DT<P>::calcGLCM( std::vector<std::vector<double> > & mat, const std::vector<int> dis, uint channel )
+{
+  assert( channel < numChannels );
+  assert( data[channel] != NULL );
+  assert( dis.size() > 2 );
+  
+  P min, max;
+  statistics( min, max, channel );
+  
+  long k = 0;
+
+  mat = std::vector<std::vector<double> > ( max+1, std::vector<double> ( max+1, 0.0 ) );
+  
+  for (int z = 0; z < zsize; z++ )
+  {
+    for (int y = 0; y < ysize; y++ )
+    {
+      for (int x = 0; x < xsize; x++)
+      {
+        if ( (x+dis[0]>=0) && (x+dis[0]<xsize)
+          && (y+dis[1]>=0) && (y+dis[1]<ysize)
+          && (z+dis[2]>=0) && (z+dis[2]<zsize) )
+        {
+          P val = get( x, y, z, channel );
+          P dval = get( x+dis[0], y+dis[1], z+dis[2], channel );
+          mat[val][dval]++;
+          mat[dval][val]++;   // symmetry
+          k += 2;
+        }
+      }
+    }
+  }
+  
+  for (int i = 0; i <= max; i++)
+  {
+    for (int j = 0; j <= max; j++)
+    {
+      mat[i][j] /= k;
+    }
+  }
+
+}
+
+template<class P>
+void MultiChannelImage3DT<P>::calcVariance( uint srcchan, uint tarchan )
+{
+  assert( srcchan < tarchan );
+  assert( tarchan < numChannels );
+  assert( data[srcchan] != NULL );
+  assert( data[tarchan] != NULL );
+  
+  uint windowsize = 3;
+  int win = (windowsize-1)/2;
+  
+  for ( int z = 0; z < zsize; z++ )
+  {
+    for ( int y = 0; y < ysize; y++ )
+    {
+      for ( int x = 0; x < xsize; x++ )
+      {
+        int meansum = 0;
+          for ( int u = -win; u <= win; u++ ) 
+          {
+            for ( int v = -win; v <= win; v++ )
+            {
+              for ( int w = -win; w <= win; w++)
+              {
+                int u_tmp = u;
+                int v_tmp = v;
+                int w_tmp = w;
+                if ( (x+u<0) || (x+u>=xsize) )
+                  u_tmp = -u_tmp;
+                if ( (y+v<0) || (y+v>=ysize) )
+                  v_tmp = -v_tmp;
+                if ( (z+w<0) || (z+w>=zsize) )
+                  w_tmp = -w_tmp;
+                meansum += get( x+u_tmp, y+v_tmp, z+w_tmp, srcchan );
+              }
+            }
+          }
+        meansum /= (windowsize*windowsize*windowsize);
+        
+        unsigned long varsum = 0;
+        for ( int u = -win; u <= win; u++ ) 
+        {
+          for ( int v = -win; v <= win; v++ )
+          {
+            for ( int w = -win; w <= win; w++)
+            {
+              int u_tmp = u;
+              int v_tmp = v;
+              int w_tmp = w;
+              if ( (x+u<0) || (x+u>=xsize) )
+                u_tmp = -u_tmp;
+              if ( (y+v<0) || (y+v>=ysize) )
+                v_tmp = -v_tmp;
+              if ( (z+w<0) || (z+w>=zsize) )
+                  w_tmp = -w_tmp;
+              
+              long sdev = (get( x+u_tmp, y+v_tmp, z+w_tmp, srcchan ) - meansum );
+              varsum += (sdev*sdev);
+            }
+          }
+        }
+        varsum /= (windowsize*windowsize+windowsize)-1;
+        set( x, y, z, varsum, tarchan );
+      }
+    }
+  }
+}
+
+template<class P>
+P MultiChannelImage3DT<P>::getIntegralValue(int ulfx, int ulfy, int ulfz, int lrbx, int lrby, int lrbz, int channel)
+{
+  ulfx = std::max(ulfx-1, -1);
+  ulfx = std::min(ulfx, xsize-1);
+  ulfy = std::max(ulfy-1, -1);
+  ulfy = std::min(ulfy, ysize-1);
+  ulfz = std::max(ulfz-1, -1);
+  ulfz = std::min(ulfz, zsize-1);
+
+  lrbx = std::max(lrbx, 0);
+  lrbx = std::min(lrbx, xsize-1);
+  lrby = std::max(lrby, 0);
+  lrby = std::min(lrby, ysize-1);
+  lrbz = std::max(lrbz, 0);
+  lrbz = std::min(lrbz, zsize-1);
+  
+  P val1, val2, val3, val4, val5, val6, val7, val8;
+  val1 = get(lrbx, lrby, lrbz, channel);
+  if( ulfz > -1 )
+    val2 = get(lrbx, lrby, ulfz, channel);
+  else
+    val2 = 0;
+  if( ulfx > -1 )
+    val3 = get(ulfx, lrby, lrbz, channel);
+  else
+    val3 = 0;
+  if( ulfx > -1 && ulfz > -1 )
+    val4 = get(ulfx, lrby, ulfz, channel);
+  else
+    val4 = 0;
+  if( ulfy > -1 )
+    val5 = get(lrbx, ulfy, lrbz, channel);
+  else
+    val5 = 0;
+  if( ulfy > -1 && ulfz > -1 )
+    val6 = get(lrbx, ulfy, ulfz, channel);
+  else
+    val6 = 0;
+  if( ulfx > -1 && ulfy > -1 )
+    val7 = get(ulfx, ulfy, lrbz, channel);
+  else
+    val7 = 0;
+  if(ulfx > -1 && ulfy > -1 && ulfz > -1)
+    val8 = get(ulfx, ulfy, ulfz, channel);
+  else
+    val8 = 0;
+  
+  P volume = abs((P)(lrbx-ulfx)*(lrby-ulfy)*(lrbz-ulfz));
+  P val = val1 - val2 - val3 + val4 - ( val5 - val6 - val7 + val8 );
+  return val/volume;
+}
+
+template<class P>
+void MultiChannelImage3DT<P>::store( std::string filename ) const
+{
+  // simple raw format
+  FILE *f = fopen( filename.c_str(), "w" );
+
+  if ( f == NULL ) {
+    fprintf( stderr, "MultiChannelImage3DT::store: error writing to %s\n", filename.c_str() );
+    exit( -1 );
+  }
+
+  fwrite( &xsize, sizeof( int ), 1, f );
+  fwrite( &ysize, sizeof( int ), 1, f );
+  fwrite( &zsize, sizeof( int ), 1, f );
+  fwrite( &numChannels, sizeof( uint ), 1, f );
+
+  for ( uint channel = 0 ; channel < numChannels ; channel++ )
+  {
+    assert( data[channel] != NULL );
+    fwrite( data[channel], sizeof( P ), xsize*ysize*zsize, f );
+  }
+
+  fclose( f );
+}
+
+template<class P>
+void MultiChannelImage3DT<P>::restore( std::string filename )
+{
+  // simple raw format
+  FILE *f = fopen( filename.c_str(), "r" );
+
+  if ( f == NULL ) {
+    fprintf( stderr, "MultiChannelImage3DT::store: error reading from %s\n", filename.c_str() );
+    exit( -1 );
+  }
+
+  fread( &xsize, sizeof( int ), 1, f );
+  fread( &ysize, sizeof( int ), 1, f );
+  fread( &zsize, sizeof( int ), 1, f );
+  fread( &numChannels, sizeof( uint ), 1, f );
+
+  if ( numChannels > 0 ) {
+    reInit( xsize, ysize, zsize, numChannels );
+
+    for ( uint channel = 0 ; channel < numChannels ; channel++ )
+    {
+      assert( data[channel] != NULL );
+      fread( data[channel], sizeof( P ), xsize*ysize*zsize, f );
+    }
+  } else {
+    freeData();
+    data = NULL;
+  }
+
+  fclose( f );
+}
+
+template<class P>
+int MultiChannelImage3DT<P>::width() const
+{
+  return xsize;
+}
+
+template<class P>
+int MultiChannelImage3DT<P>::height() const
+{
+  return ysize;
+}
+
+template<class P>
+int MultiChannelImage3DT<P>::depth() const
+{
+  return zsize;
+}
+
+template<class P>
+int MultiChannelImage3DT<P>::channels() const
+{
+  return ( int )numChannels;
+}
+
+template<class P>
+int MultiChannelImage3DT<P>::getPixelInt( int x, int y, int z, int channel ) const
+{
+  throw( "this type is not implemented\n" );
+  return -1;
+}
+
+template<class P>
+double MultiChannelImage3DT<P>::getPixelFloat( int x, int y, int z, int channel ) const
+{
+  throw( "this type is not implemented\n" );
+  return -1.0;
+}
+
+template<class P>
+void MultiChannelImage3DT<P>::setPixelInt( int x, int y, int z, int channel, int pixel )
+{
+  throw( "this type is not implemented\n" );
+}
+
+template<class P>
+void MultiChannelImage3DT<P>::setPixelFloat( int x, int y, int z, int channel, double pixel )
+{
+  throw( "this type is not implemented\n" );
+}
+
+#define SET_FUNCS_PROTO_MACRO3D(MYTYPE) \
+template<>\
+int MultiChannelImage3DT<MYTYPE>::getPixelInt(int x, int y, int z, int channel) const;\
+template<>\
+double MultiChannelImage3DT<MYTYPE>::getPixelFloat(int x, int y, int z, int channel) const;\
+template<>\
+void MultiChannelImage3DT<MYTYPE>::setPixelInt(int x, int y, int z, int channel, int pixel);\
+template<>\
+void MultiChannelImage3DT<MYTYPE>::setPixelFloat(int x, int y, int z, int channel, double pixel);
+
+SET_FUNCS_PROTO_MACRO3D( double )
+SET_FUNCS_PROTO_MACRO3D( int )
+SET_FUNCS_PROTO_MACRO3D( long int )
+SET_FUNCS_PROTO_MACRO3D( float )
+SET_FUNCS_PROTO_MACRO3D( unsigned int )
+
+} // namespace
+
+

core/image/MultiChannelImageAccess3D.h

@@ -0,0 +1,35 @@
+/*
+ * NICE-Core - efficient algebra and computer vision methods
+ *  - libimage - An image library
+ * See file License for license information.
+ */
+
+#ifndef _LIMUN_BASEMULTICHANNELIMAGE3D_H
+#define _LIMUN_BASEMULTICHANNELIMAGE3D_H
+
+namespace NICE {
+
+/**
+ * A generic, but not very efficient multi channel image access interface.
+ *
+ * @author Sven Sickert
+ */
+class MultiChannelImageAccess3D {
+public:
+  virtual int width() const = 0;
+  virtual int height() const = 0;
+	virtual int depth() const = 0;
+  virtual int channels() const = 0;
+
+  virtual int getPixelInt(int x, int y, int z, int channel) const = 0;
+  virtual double getPixelFloat(int x, int y, int z, int channel) const = 0;
+
+  virtual void setPixelInt(int x, int y, int z, int channel, int pixel) = 0;
+  virtual void setPixelFloat(int x, int y, int z, int channel, double pixel) = 0;
+
+  virtual ~MultiChannelImageAccess3D () {};
+};
+
+} // namespace
+
+#endif

core/vector/MatrixT.tcc

@@ -417,11 +417,16 @@ MatrixT<ElementType>::operator _Op##= (const MatrixT<ElementType>& v) {   \
   template<typename ElementType>
   void MatrixT<ElementType>::setRow(const uint & row, const NICE::VectorT<ElementType> & v)
   {
-    if (rows() < v.size() )
+    if (cols() < v.size() )
     {
       _THROW_EMatrix("Matrix setRow: target matrix too small.");
     }
-
+    
+    if (row >= rows() )
+    {
+      _THROW_EMatrix("Matrix setRow: row does not specify a proper row index.");
+    }
+    
     // FIXME use IPP?
     for (unsigned int j = 0; j < v.size(); j++)
     {

core/vector/tests/test_setRow.txt

@@ -0,0 +1,16 @@
+MatrixT<double> a(1,2,4.5);
+MatrixT<double> b(2,1,4.5);
+MatrixT<double> c(1,3,4.5);
+VectorT<double> v(2,1.0);
+
+a.setRow(0,v);
+//a(0,0) and a(0,1) should now be 1.0 and not 4.5
+
+a.setRow(1,v)
+//should fail, because this row is non-existent
+
+b.setRow(0,v)
+//should fail, because vector is too long 
+
+c.setRow(0,v)
+//should give a warning, because vector is shorter than one complete row