merge branch 'semseg3d' into 'master'

libdepend.inc

@@ -1 +1,2 @@
+$(call PKG_DEPEND_INT,segmentation)
 $(call PKG_DEPEND_INT,core)

progs/libdepend.inc

@@ -0,0 +1,2 @@
+$(call PKD_DEPEND_INT,segmentation)
+$(call PKG_DEPEND_INT,core)

progs/testSemanticSegmentation3D.cpp

@@ -0,0 +1,387 @@
+// Beispielhafter Aufruf: BUILD_x86_64/progs/testSemanticSegmentation -config <CONFIGFILE>
+
+/**
+* @file testSemanticSegmentation.cpp
+* @brief test semantic segmentation routines for 3d images and 2d images
+* @author Erik Rodner, Björn Fröhlich, Sven Sickert
+* @date 03/20/2008
+*/
+
+#ifdef NICE_USELIB_OPENMP
+#include <omp.h>
+#endif
+
+#include "core/basics/Config.h"
+#include "core/basics/StringTools.h"
+#include <vislearning/baselib/ICETools.h>
+
+#include <core/image/MultiChannelImage3DT.h>
+#include <semseg3d/semseg/SemSegContextTree3D.h>
+
+#include <core/basics/ResourceStatistics.h>
+#include <core/image/Morph.h>
+
+#include <fstream>
+#include <vector>
+
+#undef DEBUG
+
+using namespace OBJREC;
+
+using namespace NICE;
+
+using namespace std;
+
+void segmentToOverlay ( const NICE::Image *orig, const NICE::ColorImage & segment,
+                        NICE::ColorImage & result )
+{
+  int xsize = orig->width();
+  int ysize = orig->height();
+  
+  result.resize( xsize, ysize );
+  vector< NICE::MatrixT<double> > channelMat;
+  
+  double alpha = .5;
+  
+  for (int c = 0; c < 3; c++)
+  {
+    NICE::MatrixT<double> chan ( xsize, ysize );
+    channelMat.push_back( chan );
+  }
+
+  for (int y = 0; y < ysize; y++)
+  {
+    for (int x = 0; x < xsize; x++)
+    {
+      uchar val = orig->getPixelQuick(x,y);
+      for (int c = 0; c < 3; c++)
+        channelMat[c](x,y) = (double)val + alpha*(double)segment.getPixel( x, y, c );
+    }
+  }
+  
+  for (int c = 0; c < 3; c++)
+  {
+    channelMat[c] /= channelMat[c].Max();
+    channelMat[c] *= 255;
+  }
+  
+  for (int y = 0; y < ysize; y++)
+  {
+    for (int x = 0; x < xsize; x++)
+    {
+      for (int c = 0; c < 3; c++)
+      {
+        int val = channelMat[c](x,y);
+        result.setPixel( x, y, c, (uchar)val);
+      }
+    }
+  }
+}
+
+void updateMatrix ( const NICE::Image & img, const NICE::Image & gt,
+                    NICE::Matrix & M, const set<int> & forbidden_classes )
+{
+  double subsamplex = gt.width() / ( double ) img.width();
+  double subsampley = gt.height() / ( double ) img.height();
+
+  for ( int y = 0 ; y < gt.height() ; y++ )
+    for ( int x = 0 ; x < gt.width() ; x++ )
+    {
+      int xx = ( int ) ( x / subsamplex );
+      int yy = ( int ) ( y / subsampley );
+
+      if ( xx < 0 ) xx = 0;
+
+      if ( yy < 0 ) yy = 0;
+
+      if ( xx > img.width() - 1 ) xx = img.width() - 1;
+
+      if ( yy > img.height() - 1 ) yy = img.height() - 1;
+
+      int cimg = img.getPixel ( xx, yy );
+
+      int gimg = gt.getPixel ( x, y );
+
+      if ( forbidden_classes.find ( gimg ) == forbidden_classes.end() )
+      {
+        M ( gimg, cimg ) ++;
+      }
+    }
+}
+
+/**
+ test semantic segmentation routines
+*/
+int main ( int argc, char **argv )
+{
+  std::set_terminate ( __gnu_cxx::__verbose_terminate_handler );
+
+  Config conf ( argc, argv );
+
+  ResourceStatistics rs;
+
+  /*-------------I/O CONFIGURATION-------------*/
+  bool postProcessing = conf.gB( "main", "post_process", false);
+  bool run_3Dseg = conf.gB( "SSContextTree", "run_3dseg", false);
+  bool show_result = conf.gB ( "debug", "show_results", false );
+  bool write_results = conf.gB ( "debug", "write_results", false );
+  string output_type = conf.gS ( "debug", "output_type", "ppm" );
+  string output_postfix = conf.gS ( "debug", "output_postfix", "" );
+  string resultdir = conf.gS ( "debug", "resultdir", "." );
+  /*-------------------------------------------*/
+
+#ifdef DEBUG
+  cerr << "Writing Results to " << resultdir << endl;
+#endif
+
+
+  MultiDataset md ( &conf );
+
+  const ClassNames & classNames = md.getClassNames ( "train" );
+
+  // initialize semantic segmentation method
+  SemanticSegmentation *semseg = NULL;
+  semseg = new SemSegContextTree3D ( &conf, &md );
+  
+  // train semantic segmentation method
+  cout << "\nTRAINING" << endl;
+  cout << "########\n" << endl;
+  semseg->train( &md );
+  
+  const LabeledSet *testFiles = md["test"];
+
+  set<int> forbidden_classes;
+  std::string forbidden_classes_s = conf.gS ( "analysis", "forbidden_classes", "" );
+
+  classNames.getSelection ( forbidden_classes_s, forbidden_classes );
+
+//   ProgressBar pb ( "Semantic Segmentation Analysis" );
+//   pb.show();
+
+  vector< int > zsizeVec;
+  semseg->getDepthVector ( testFiles, zsizeVec, run_3Dseg );
+
+  int depthCount = 0, idx = 0;
+  vector< string > filelist;
+  NICE::MultiChannelImageT<double> segresult;
+  NICE::MultiChannelImageT<double> gt;
+  std::vector< NICE::Matrix > M_vec;
+
+  cout << "\nCLASSIFICATION" << endl;
+  cout << "##############\n" << endl;
+  for (LabeledSet::const_iterator it = testFiles->begin(); it != testFiles->end(); it++)
+  {
+    for (std::vector<ImageInfo *>::const_iterator jt = it->second.begin();
+         jt != it->second.end(); jt++)
+    {
+      ImageInfo & info = *(*jt);
+      std::string file = info.img();
+      filelist.push_back ( file );
+      depthCount++;
+
+      NICE::Image lm;
+      NICE::Image lm_gt;
+      if ( info.hasLocalizationInfo() )
+      {
+        const LocalizationResult *l_gt = info.localization();
+
+        lm.resize ( l_gt->xsize, l_gt->ysize );
+        lm.set ( 0 );
+
+        lm_gt.resize ( l_gt->xsize, l_gt->ysize );
+        lm_gt.set ( 0 );
+
+        l_gt->calcLabeledImage ( lm, classNames.getBackgroundClass() );
+#ifdef DEBUG
+        cout << "testSemanticSegmentation3D: Generating Labeled NICE::Image (Ground-Truth)" << endl;
+#endif
+        l_gt->calcLabeledImage ( lm_gt, classNames.getBackgroundClass() );
+      }
+      segresult.addChannel ( lm );
+      gt.addChannel ( lm_gt );
+
+      int depthBoundary = 0;
+      if ( run_3Dseg )
+      {
+        depthBoundary = zsizeVec[idx];
+      }
+
+      if ( depthCount < depthBoundary ) continue;
+
+      NICE::MultiChannelImage3DT<double> probabilities;
+      NICE::MultiChannelImage3DT<double> imgData;
+
+      semseg->make3DImage ( filelist, imgData );
+      semseg->classify ( imgData, segresult, probabilities, filelist );
+
+      // save to file
+      for ( int z = 0; z < segresult.channels(); z++ )
+      {
+        std::string fname = StringTools::baseName ( filelist[z], false );
+
+        if ( show_result || write_results )
+        {
+          NICE::ColorImage orig ( filelist[z] );
+          NICE::ColorImage rgb;
+          NICE::ColorImage rgb_gt;
+          NICE::ColorImage ov_rgb;
+          NICE::ColorImage ov_rgb_gt;
+
+          for ( int y = 0 ; y < segresult.height(); y++ )
+          {
+            for ( int x = 0 ; x < segresult.width(); x++ )
+            {
+              lm.setPixel ( x, y, segresult.get ( x, y, ( uint ) z ) );
+              if ( run_3Dseg )
+                lm_gt.setPixel ( x, y, gt.get ( x, y, ( uint ) z ) );
+            }
+          }
+
+          // confusion matrix
+          NICE::Matrix M ( classNames.getMaxClassno() + 1, classNames.getMaxClassno() + 1 );
+          M.set ( 0 );
+          updateMatrix ( lm, lm_gt, M, forbidden_classes );
+          M_vec.push_back ( M );
+
+          classNames.labelToRGB ( lm, rgb );
+          classNames.labelToRGB ( lm_gt, rgb_gt );
+
+          if (postProcessing)
+          {
+            // median filter
+            for (int r = 0; r < 3; r++)
+            {
+              NICE::Image postIm(rgb.width(), rgb.height());
+              NICE::median(*(rgb.getChannel(r)), &postIm, 1);
+              for (int y = 0; y < rgb.height(); y++)
+                for (int x = 0; x < rgb.width(); x++)
+                  rgb.setPixel(x,y,r, postIm.getPixelQuick(x,y));
+            }
+          }
+
+          segmentToOverlay ( orig.getChannel(1), rgb, ov_rgb );
+          segmentToOverlay ( orig.getChannel(1), rgb_gt, ov_rgb_gt );
+          
+          if ( write_results )
+          {
+            std::stringstream out;       
+            if ( output_postfix.size() > 0 )
+              out << resultdir << "/" << fname << output_postfix;
+            else
+              out << resultdir << "/" << fname;
+#ifdef DEBUG
+            cout << "Writing to file " << out.str() << "_*." << output_type << endl;
+#endif
+            orig.write ( out.str() + "_orig." + output_type );
+            rgb.write ( out.str() + "_result." + output_type );
+            rgb_gt.write ( out.str() + "_groundtruth." + output_type );
+            ov_rgb.write ( out.str() + "_overlay_res." + output_type );
+            ov_rgb_gt.write ( out.str() + "_overlay_gt." + output_type );
+          }
+        }
+      }
+
+      // prepare for new 3d image
+      filelist.clear();
+      segresult.reInit(0,0,0);
+      gt.reInit(0,0,0);
+      depthCount = 0;
+      idx++;
+//       pb.update ( testFiles->count() );
+    }
+  }
+  
+  segresult.freeData();
+//   pb.hide();
+
+  cout << "\nSTATISTICS" << endl;
+  cout << "##########\n" << endl;
+  long maxMemory;
+  double userCPUTime, sysCPUTime;
+  rs.getStatistics ( maxMemory, userCPUTime, sysCPUTime );
+  cout << "Memory (max):    " << maxMemory << " KB" << endl;
+  cout << "CPU Time (user): " << userCPUTime << " seconds" << endl;
+  cout << "CPU Time (sys):  " << sysCPUTime << " seconds" << endl;
+
+  double overall = 0.0;
+  double sumall = 0.0;
+
+  NICE::Matrix M ( classNames.getMaxClassno() + 1, classNames.getMaxClassno() + 1 );
+  M.set ( 0 );
+  for ( int s = 0; s < ( int ) M_vec.size(); s++ )
+  {
+    NICE::Matrix M_tmp = M_vec[s];
+    for ( int r = 0; r < ( int ) M_tmp.rows(); r++ )
+    {
+      for ( int c = 0; c < ( int ) M_tmp.cols(); c++ )
+      {
+        if ( r == c )
+          overall += M_tmp ( r, c );
+
+        sumall += M_tmp ( r, c );
+        M ( r, c ) += M_tmp ( r, c );
+      }
+    }
+  }
+  overall /= sumall;
+
+  // normalizing M using rows
+  for ( int r = 0 ; r < ( int ) M.rows() ; r++ )
+  {
+    double sum = 0.0;
+
+    for ( int c = 0 ; c < ( int ) M.cols() ; c++ )
+      sum += M ( r, c );
+
+    if ( fabs ( sum ) > 1e-4 )
+      for ( int c = 0 ; c < ( int ) M.cols() ; c++ )
+        M ( r, c ) /= sum;
+  }
+
+  double avg_perf = 0.0;
+  int classes_trained = 0;
+
+  for ( int r = 0 ; r < ( int ) M.rows() ; r++ )
+  {
+    if ( ( classNames.existsClassno ( r ) ) && ( forbidden_classes.find ( r ) == forbidden_classes.end() ) )
+    {
+      avg_perf += M ( r, r );
+      double lsum = 0.0;
+      for ( int r2 = 0; r2 < ( int ) M.rows(); r2++ )
+      {
+        lsum += M ( r,r2 );
+      }
+      if ( lsum != 0.0 )
+      {
+        classes_trained++;
+      }
+    }
+  }
+
+  // print/save results of evaluation
+  cout << "\nPERFORMANCE" << endl;
+  cout << "###########\n" << endl;
+  ofstream fout ( ( resultdir + "/res.txt" ).c_str(), ios::out );
+  fout << "Overall Recognition Rate: " << overall << endl;
+  fout << "Average Recognition Rate: " << avg_perf / ( classes_trained ) << endl;
+  fout << "Lower Bound: " << 1.0  / classes_trained << endl;
+  cout << "Overall Recogntion Rate: " << overall << endl;
+  cout << "Average Recogntion Rate: " << avg_perf / ( classes_trained ) << endl;
+  cout << "Lower Bound: " << 1.0  / classes_trained << endl;
+
+  cout <<"\nClasses:" << endl;
+  for ( int r = 0 ; r < ( int ) M.rows() ; r++ )
+  {
+    if ( ( classNames.existsClassno ( r ) ) && ( forbidden_classes.find ( r ) == forbidden_classes.end() ) )
+    {
+      std::string classname = classNames.text ( r );
+      fout << classname.c_str() << ": " << M ( r, r ) << endl;
+      cout << classname.c_str() << ": " << M ( r, r ) << endl;
+    }
+  }
+  fout.close();
+
+  delete semseg;
+
+  return 0;
+}

semseg/SemSegContextTree3D.cpp

@@ -0,0 +1,2279 @@
+#include "SemSegContextTree3D.h"
+#include "vislearning/baselib/Globals.h"
+#include "core/basics/StringTools.h"
+
+#include "core/imagedisplay/ImageDisplay.h"
+
+#include "vislearning/cbaselib/CachedExample.h"
+#include "vislearning/cbaselib/PascalResults.h"
+#include "vislearning/baselib/cc.h"
+#include "segmentation/RSMeanShift.h"
+#include "segmentation/RSGraphBased.h"
+#include "segmentation/RSSlic.h"
+#include "core/basics/numerictools.h"
+#include "core/basics/StringTools.h"
+#include "core/basics/FileName.h"
+#include "vislearning/baselib/ICETools.h"
+
+#include "core/basics/Timer.h"
+#include "core/basics/vectorio.h"
+#include "core/basics/quadruplet.h"
+#include <core/image/Filter.h>
+#include "core/image/FilterT.h"
+#include <core/image/Morph.h>
+
+#include <omp.h>
+#include <iostream>
+
+#define VERBOSE
+#undef DEBUG
+#undef VISUALIZE
+#undef WRITEREGIONS
+
+using namespace OBJREC;
+using namespace std;
+using namespace NICE;
+
+
+//###################### CONSTRUCTORS #########################//
+
+
+SemSegContextTree3D::SemSegContextTree3D () : SemanticSegmentation ()
+{
+  this->lfcw                = NULL;
+  this->firstiteration      = true;
+  this->run3Dseg            = false;
+  this->maxSamples          = 2000;
+  this->minFeats            = 50;
+  this->maxDepth            = 10;
+  this->windowSize          = 15;
+  this->contextMultiplier   = 3;
+  this->featsPerSplit       = 200;
+  this->useShannonEntropy   = true;
+  this->nbTrees             = 10;
+  this->randomTests         = 10;
+  this->useAltTristimulus   = false;
+  this->useGradient         = true;
+  this->useWeijer           = false;
+  this->useAdditionalLayer  = false;
+  this->useHoiemFeatures    = false;
+  this->useCategorization   = false;
+  this->cndir               = "";
+  this->fasthik             = NULL;
+  this->saveLoadData        = false;
+  this->fileLocation        = "tmp.txt";
+  this->pixelWiseLabeling   = true;
+  this->segmentation        = NULL;
+  this->useFeat0            = true;
+  this->useFeat1            = false;
+  this->useFeat2            = true;
+  this->useFeat3            = true;
+  this->useFeat4            = false;
+  this->useFeat5            = false;
+}
+
+
+SemSegContextTree3D::SemSegContextTree3D (
+    const Config *conf,
+    const MultiDataset *md )
+    : SemanticSegmentation ( conf, & ( md->getClassNames ( "train" ) ) )
+{
+  this->conf = conf;
+
+  string section = "SSContextTree";
+  string featsec = "Features";
+
+  this->lfcw                = NULL;
+  this->firstiteration      = true;
+  this->run3Dseg            = conf->gB ( section, "run_3dseg", false );
+  this->maxSamples          = conf->gI ( section, "max_samples", 2000 );
+  this->minFeats            = conf->gI ( section, "min_feats", 50 );
+  this->maxDepth            = conf->gI ( section, "max_depth", 10 );
+  this->windowSize          = conf->gI ( section, "window_size", 15 );
+  this->contextMultiplier   = conf->gI ( section, "context_multiplier", 3 );
+  this->featsPerSplit       = conf->gI ( section, "feats_per_split", 200 );
+  this->useShannonEntropy   = conf->gB ( section, "use_shannon_entropy", true );
+  this->nbTrees             = conf->gI ( section, "amount_trees", 10 );
+  this->randomTests         = conf->gI ( section, "random_tests", 10 );
+
+  this->useAltTristimulus   = conf->gB ( featsec, "use_alt_trist", false );
+  this->useGradient         = conf->gB ( featsec, "use_gradient", true );
+  this->useWeijer           = conf->gB ( featsec, "use_weijer", true );
+  this->useAdditionalLayer  = conf->gB ( featsec, "use_additional_layer", false );
+  this->useHoiemFeatures    = conf->gB ( featsec, "use_hoiem_features", false );
+
+  this->useCategorization   = conf->gB ( section, "use_categorization", false );
+  this->cndir               = conf->gS ( "SSContextTree", "cndir", "" );
+
+  this->saveLoadData        = conf->gB ( "debug", "save_load_data", false );
+  this->fileLocation        = conf->gS ( "debug", "datafile", "tmp.txt" );
+
+  this->pixelWiseLabeling   = conf->gB ( section, "pixelWiseLabeling", false );
+
+  if ( useCategorization && cndir == "" )
+    this->fasthik = new GPHIKClassifierNICE ( conf );
+  else
+    this->fasthik = NULL;
+
+  if ( useWeijer )
+    this->lfcw    = new LocalFeatureColorWeijer ( conf );
+
+  this->classnames = md->getClassNames ( "train" );
+
+  // feature types
+  this->useFeat0 = conf->gB ( section, "use_feat_0", true);     // pixel pair features
+  this->useFeat1 = conf->gB ( section, "use_feat_1", false);    // region feature
+  this->useFeat2 = conf->gB ( section, "use_feat_2", true);     // integral features
+  this->useFeat3 = conf->gB ( section, "use_feat_3", true);     // integral contex features
+  this->useFeat4 = conf->gB ( section, "use_feat_4", false);    // pixel pair context features
+  this->useFeat5 = conf->gB ( section, "use_feat_5", false);    // ray features
+
+  string segmentationtype = conf->gS ( section, "segmentation_type", "slic" );
+  if ( segmentationtype == "meanshift" )
+    this->segmentation = new RSMeanShift ( conf );
+  else if ( segmentationtype == "felzenszwalb" )
+    this->segmentation = new RSGraphBased ( conf );
+  else if ( segmentationtype == "slic" )
+    this->segmentation = new RSSlic ( conf );
+  else if ( segmentationtype == "none" )
+  {
+    this->segmentation = NULL;
+    this->pixelWiseLabeling = true;
+    this->useFeat1 = false;
+  }
+  else
+    throw ( "no valid segmenation_type\n please choose between none, meanshift, slic and felzenszwalb\n" );
+
+  if ( !useGradient || imagetype == IMAGETYPE_RGB)
+    this->useFeat5 = false;
+
+  if ( useFeat0 )
+    this->featTypes.push_back(0);
+  if ( useFeat1 )
+    this->featTypes.push_back(1);
+  if ( useFeat2 )
+    this->featTypes.push_back(2);
+  if ( useFeat3 )
+    this->featTypes.push_back(3);
+  if ( useFeat4 )
+    this->featTypes.push_back(4);
+  if ( useFeat5 )
+    this->featTypes.push_back(5);
+
+  this->initOperations();
+}
+
+//###################### DESTRUCTORS ##########################//
+
+SemSegContextTree3D::~SemSegContextTree3D()
+{
+}
+
+
+//#################### MEMBER FUNCTIONS #######################//
+
+void SemSegContextTree3D::initOperations()
+{
+  string featsec = "Features";
+
+  // operation prototypes
+  vector<Operation3D*> tops0, tops1, tops2, tops3, tops4, tops5;
+
+  if ( conf->gB ( featsec, "int", true ) )
+  {
+    tops2.push_back ( new IntegralOps() );
+    Operation3D* o = new IntegralOps();
+    o->setContext(true);
+    tops3.push_back ( o );
+  }
+  if ( conf->gB ( featsec, "bi_int_cent", true ) )
+  {
+    tops2.push_back ( new BiIntegralCenteredOps() );
+    Operation3D* o = new BiIntegralCenteredOps();
+    o->setContext(true);
+    tops3.push_back ( o );
+  }
+  if ( conf->gB ( featsec, "int_cent", true ) )
+  {
+    tops2.push_back ( new IntegralCenteredOps() );
+    Operation3D* o = new IntegralCenteredOps();
+    o->setContext(true);
+    tops3.push_back ( o );
+  }
+  if ( conf->gB ( featsec, "haar_horz", true ) )
+  {
+    tops2.push_back ( new HaarHorizontal() );
+    Operation3D* o = new HaarHorizontal();
+    o->setContext(true);
+    tops3.push_back ( o );
+  }
+  if ( conf->gB ( featsec, "haar_vert", true ) )
+  {
+    tops2.push_back ( new HaarVertical );
+    Operation3D* o = new HaarVertical();
+    o->setContext(true);
+    tops3.push_back ( o );
+  }
+  if ( conf->gB ( featsec, "haar_stack", true ) )
+  {
+    tops2.push_back ( new HaarStacked() );
+    Operation3D* o = new HaarStacked();
+    o->setContext(true);
+    tops3.push_back ( o );
+  }
+  if ( conf->gB ( featsec, "haar_diagxy", true ) )
+  {
+    tops2.push_back ( new HaarDiagXY() );
+    Operation3D* o = new HaarDiagXY();
+    o->setContext(true);
+    tops3.push_back ( o );
+  }
+  if ( conf->gB ( featsec, "haar_diagxz", true ) )
+  {
+    tops2.push_back ( new HaarDiagXZ() );
+    Operation3D* o = new HaarDiagXZ();
+    o->setContext(true);
+    tops3.push_back ( o );
+  }
+  if ( conf->gB ( featsec, "haar_diagyz", true ) )
+  {
+    tops2.push_back ( new HaarDiagYZ() );
+    Operation3D* o = new HaarDiagYZ();
+    o->setContext(true);
+    tops3.push_back ( o );
+  }
+  if ( conf->gB ( featsec, "haar3_horz", true ) )
+  {
+    tops2.push_back ( new Haar3Horiz() );
+    Operation3D* o = new Haar3Horiz();
+    o->setContext(true);
+    tops3.push_back ( o );
+  }
+  if ( conf->gB ( featsec, "haar3_vert", true ) )
+  {
+    tops2.push_back ( new Haar3Vert() );
+    Operation3D* o = new Haar3Vert();
+    o->setContext(true);
+    tops3.push_back ( o );
+  }
+  if ( conf->gB ( featsec, "haar3_stack", true ) )
+  {
+    tops2.push_back ( new Haar3Stack() );
+    Operation3D* o = new Haar3Stack();
+    o->setContext(true);
+    tops3.push_back ( o );
+  }
+
+  if ( conf->gB ( featsec, "minus", true ) )
+  {
+    tops0.push_back ( new Minus() );
+    Operation3D* o = new Minus();
+    o->setContext(true);
+    tops4.push_back ( o );
+  }
+  if ( conf->gB ( featsec, "minus_abs", true ) )
+  {
+    tops0.push_back ( new MinusAbs() );
+    Operation3D* o = new MinusAbs();
+    o->setContext(true);
+    tops4.push_back ( o );
+  }
+  if ( conf->gB ( featsec, "addition", true ) )
+  {
+    tops0.push_back ( new Addition() );
+    Operation3D* o = new Addition();
+    o->setContext(true);
+    tops4.push_back ( o );
+  }
+  if ( conf->gB ( featsec, "only1", true ) )
+  {
+    tops0.push_back ( new Only1() );
+    Operation3D* o = new Only1();
+    o->setContext(true);
+    tops4.push_back ( o );
+  }
+  if ( conf->gB ( featsec, "rel_x", true ) )
+    tops0.push_back ( new RelativeXPosition() );
+  if ( conf->gB ( featsec, "rel_y", true ) )
+    tops0.push_back ( new RelativeYPosition() );
+  if ( conf->gB ( featsec, "rel_z", true ) )
+    tops0.push_back ( new RelativeZPosition() );
+
+  if ( conf->gB ( featsec, "ray_diff", false ) )
+    tops5.push_back ( new RayDiff() );
+  if ( conf->gB ( featsec, "ray_dist", false ) )
+    tops5.push_back ( new RayDist() );
+  if ( conf->gB ( featsec, "ray_orient", false ) )
+    tops5.push_back ( new RayOrient() );
+  if ( conf->gB ( featsec, "ray_norm", false ) )
+    tops5.push_back ( new RayNorm() );
+
+  this->ops.push_back ( tops0 );
+  this->ops.push_back ( tops1 );
+  this->ops.push_back ( tops2 );
+  this->ops.push_back ( tops3 );
+  this->ops.push_back ( tops4 );
+  this->ops.push_back ( tops5 );
+}
+
+double SemSegContextTree3D::getBestSplit (
+    std::vector<NICE::MultiChannelImage3DT<double> > &feats,
+    std::vector<NICE::MultiChannelImage3DT<unsigned short int> > &nodeIndices,
+    const std::vector<NICE::MultiChannelImageT<int> > &labels,
+    int node,
+    Operation3D *&splitop,
+    double &splitval,
+    const int &tree,
+    vector<vector<vector<double> > > &regionProbs )
+{
+  Timer t;
+  t.start();
+  int imgCount = 0;
+
+  try
+  {
+    imgCount = ( int ) feats.size();
+  }
+  catch ( Exception )
+  {
+    cerr << "no features computed?" << endl;
+  }
+
+  double bestig = -numeric_limits< double >::max();
+
+  splitop = NULL;
+  splitval = -1.0;
+
+  vector<quadruplet<int,int,int,int> > selFeats;
+  map<int, int> e;
+  int featcounter = forest[tree][node].featcounter;
+
+  if ( featcounter < minFeats )
+  {
+    return 0.0;
+  }
+
+  vector<double> fraction ( a.size(), 0.0 );
+
+  for ( uint i = 0; i < fraction.size(); i++ )
+  {
+    if ( forbidden_classes.find ( labelmapback[i] ) != forbidden_classes.end() )
+      fraction[i] = 0;
+    else
+      fraction[i] = ( ( double ) maxSamples ) / ( ( double ) featcounter * a[i] * a.size() );
+  }
+
+  featcounter = 0;
+
+  for ( int iCounter = 0; iCounter < imgCount; iCounter++ )
+  {
+    int xsize = ( int ) nodeIndices[iCounter].width();
+    int ysize = ( int ) nodeIndices[iCounter].height();
+    int zsize = ( int ) nodeIndices[iCounter].depth();
+
+    for ( int x = 0; x < xsize; x++ )
+      for ( int y = 0; y < ysize; y++ )
+        for ( int z = 0; z < zsize; z++ )
+        {
+          if ( nodeIndices[iCounter].get ( x, y, z, tree ) == node )
+          {
+            int cn = labels[iCounter].get ( x, y, ( uint ) z );
+            double randD = ( double ) rand() / ( double ) RAND_MAX;
+
+            if ( labelmap.find ( cn ) == labelmap.end() )
+              continue;
+
+            if ( randD < fraction[labelmap[cn]] )
+            {
+              quadruplet<int,int,int,int> quad( iCounter, x, y, z );
+              featcounter++;
+              selFeats.push_back ( quad );
+              e[cn]++;
+            }
+          }
+        }
+  }
+
+  // global entropy
+  double globent = 0.0;
+  for ( map<int, int>::iterator mapit = e.begin() ; mapit != e.end(); mapit++ )
+  {
+    double p = ( double ) ( *mapit ).second / ( double ) featcounter;
+    globent += p * log2 ( p );
+  }
+  globent = -globent;
+
+  if ( globent < 0.5 )
+    return 0.0;
+
+  // pointers to all randomly chosen features
+  std::vector<Operation3D*> featsel;
+
+  for ( int i = 0; i < featsPerSplit; i++ )
+  {
+    int x1, x2, y1, y2, z1, z2, ft;
+
+    do
+    {
+      ft = ( int ) ( rand() % featTypes.size() );
+      ft = featTypes[ft];
+    }
+    while ( channelsPerType[ft].size() == 0 );
+
+    int tmpws = windowSize;
+
+    if ( ft == 2 || ft == 4 )
+    {
+      //use larger window size for context features
+      tmpws *= contextMultiplier;
+    }
+
+    // use region feature only with reasonable pre-segmentation
+//    if ( ft == 1 && depth < 8 )
+//    {
+//      ft = 0;
+//    }
+
+    /* random window positions */
+    double z_ratio = conf->gB ( "SSContextTree", "z_ratio", 1.0 );
+    int tmp_z =  ( int ) floor( (tmpws * z_ratio) + 0.5 );
+    double y_ratio = conf->gB ( "SSContextTree", "y_ratio", 1.0 );
+    int tmp_y =  ( int ) floor( (tmpws * y_ratio) + 0.5 );
+    x1 = ( int ) ( rand() % tmpws ) - tmpws / 2 ;
+    x2 = ( int ) ( rand() % tmpws ) - tmpws / 2 ;
+    y1 = ( int ) ( rand() % tmp_y ) - tmp_y / 2 ;
+    y2 = ( int ) ( rand() % tmp_y ) - tmp_y / 2 ;
+    z1 = ( int ) ( rand() % tmp_z ) - tmp_z / 2 ;
+    z2 = ( int ) ( rand() % tmp_z ) - tmp_z / 2 ;
+
+    // use z1/z2 as directions (angles) in ray features
+    if ( ft == 5 )
+    {
+      z1 = ( int ) ( rand() % 8 );
+      z2 = ( int ) ( rand() % 8 );
+    }
+
+//    if (conf->gB ( "SSContextTree", "z_negative_only", false ))
+//    {
+//      z1 = -abs(z1);
+//      z2 = -abs(z2);
+//    }
+
+    /* random feature maps (channels) */
+    int f1, f2;
+    f1 = ( int ) ( rand() % channelsPerType[ft].size() );
+    if ( (rand() % 2) == 0 )
+      f2 = ( int ) ( rand() % channelsPerType[ft].size() );
+    else
+      f2 = f1;
+    f1 = channelsPerType[ft][f1];
+    f2 = channelsPerType[ft][f2];
+
+    if ( ft == 1 )
+    {
+      int classes = ( int ) regionProbs[0][0].size();
+      f2 = ( int ) ( rand() % classes );
+    }
+
+    /* random extraction method (operation) */
+    int o = ( int ) ( rand() % ops[ft].size() );
+
+    Operation3D *op = ops[ft][o]->clone();
+    op->set ( x1, y1, z1, x2, y2, z2, f1, f2, ft );
+
+    if ( ft == 3 || ft == 4 )
+      op->setContext ( true );
+    else
+      op->setContext ( false );
+
+    featsel.push_back ( op );
+  }
+
+  // do actual split tests
+  for ( int f = 0; f < featsPerSplit; f++ )
+  {
+    double l_bestig = -numeric_limits< double >::max();
+    double l_splitval = -1.0;
+    vector<double> vals;
+
+    double maxval = -numeric_limits<double>::max();
+    double minval = numeric_limits<double>::max();
+    int counter = 0;
+    for ( vector<quadruplet<int,int,int,int> >::const_iterator it = selFeats.begin();
+         it != selFeats.end(); it++ )
+    {
+      Features feat;
+      feat.feats = &feats[ ( *it ).first ];
+      feat.rProbs = &regionProbs[ ( *it ).first ];
+
+      assert ( forest.size() > ( uint ) tree );
+      assert ( forest[tree][0].dist.size() > 0 );
+
+      double val = 0.0;
+      val = featsel[f]->getVal ( feat, ( *it ).second, ( *it ).third, ( *it ).fourth );
+      if ( !isfinite ( val ) )
+      {
+#ifdef DEBUG
+        cerr << "feat " << feat.feats->width() << " " << feat.feats->height() << " " << feat.feats->depth() << endl;
+        cerr << "non finite value " << val << " for " << featsel[f]->writeInfos() <<  endl << (*it).second << " " <<  (*it).third << " " << (*it).fourth << endl;
+#endif
+        val = 0.0;
+      }
+      vals.push_back ( val );
+      maxval = std::max ( val, maxval );
+      minval = std::min ( val, minval );
+    }
+
+    if ( minval == maxval )
+      continue;
+
+    // split values
+    for ( int run = 0 ; run < randomTests; run++ )
+    {
+      // choose threshold randomly
+      double sval = 0.0;
+      sval = ( (double) rand() / (double) RAND_MAX*(maxval-minval) ) + minval;
+
+      map<int, int> eL, eR;
+      int counterL = 0, counterR = 0;
+      counter = 0;
+
+      for ( vector<quadruplet<int,int,int,int> >::const_iterator it2 = selFeats.begin();
+            it2 != selFeats.end(); it2++, counter++ )
+      {
+        int cn = labels[ ( *it2 ).first ].get ( ( *it2 ).second, ( *it2 ).third, ( *it2 ).fourth );
+        //cout << "vals[counter2] " << vals[counter2] << " val: " <<  val << endl;
+
+        if ( vals[counter] < sval )
+        {
+          //left entropie:
+          eL[cn] = eL[cn] + 1;
+          counterL++;
+        }
+        else
+        {
+          //right entropie:
+          eR[cn] = eR[cn] + 1;
+          counterR++;
+        }
+      }
+
+      double leftent = 0.0;
+      for ( map<int, int>::iterator mapit = eL.begin() ; mapit != eL.end(); mapit++ )
+      {
+        double p = ( double ) ( *mapit ).second / ( double ) counterL;
+        leftent -= p * log2 ( p );
+      }
+
+      double rightent = 0.0;
+      for ( map<int, int>::iterator mapit = eR.begin() ; mapit != eR.end(); mapit++ )
+      {
+        double p = ( double ) ( *mapit ).second / ( double ) counterR;
+        rightent -= p * log2 ( p );
+      }
+
+      //cout << "rightent: " << rightent << " leftent: " << leftent << endl;
+
+      double pl = ( double ) counterL / ( double ) ( counterL + counterR );
+
+      //information gain
+      double ig = globent - ( 1.0 - pl ) * rightent - pl * leftent;
+
+      //double ig = globent - rightent - leftent;
+
+      if ( useShannonEntropy )
+      {
+        double esplit = - ( pl * log ( pl ) + ( 1 - pl ) * log ( 1 - pl ) );
+        ig = 2 * ig / ( globent + esplit );
+      }
+
+      if ( ig > l_bestig )
+      {
+        l_bestig = ig;
+        l_splitval = sval;
+      }
+    }
+
+    if ( l_bestig > bestig )
+    {
+      bestig = l_bestig;
+      splitop = featsel[f];
+      splitval = l_splitval;
+    }
+  }
+
+#ifdef DEBUG
+  cout << "globent: " << globent <<  " bestig " << bestig << " splitval: " << splitval << endl;
+#endif
+  return bestig;
+}
+
+inline double SemSegContextTree3D::getMeanProb (
+    const int &x,
+    const int &y,
+    const int &z,
+    const int &channel,
+    const MultiChannelImage3DT<unsigned short int> &nodeIndices )
+{
+  double val = 0.0;
+
+  for ( int tree = 0; tree < nbTrees; tree++ )
+  {
+    val += forest[tree][nodeIndices.get ( x,y,z,tree ) ].dist[channel];
+  }
+
+  return val / ( double ) nbTrees;
+}
+
+void SemSegContextTree3D::computeRayFeatImage (
+    NICE::MultiChannelImage3DT<double> &feats,
+    int firstChannel )
+{
+  int xsize = feats.width();
+  int ysize = feats.height();
+  int zsize = feats.depth();
+
+  const int amountDirs = 8;
+
+  // compute ray feature maps from canny image
+  for ( int z = 0; z < zsize; z++)
+  {
+    // canny image from raw channel
+    NICE::Image med (xsize,ysize);
+    NICE::median ( feats.getChannel( z, 0 ), &med, 2);
+    NICE::Image* can = NICE::canny( med, 5, 25);
+
+    for ( int dir = 0; dir < amountDirs; dir++)
+    {
+      NICE::Matrix dist(xsize,ysize,0);
+      NICE::Matrix norm(xsize,ysize,0);
+      NICE::Matrix orient(xsize,ysize,0);
+
+      for (int y = 0; y < ysize; y++)
+        for ( int x = 0; x < xsize; x++)
+        {
+          int xo = 0, yo = 0; // offsets
+          int theta = 0;
+
+          switch (dir)
+          {
+            case 0: theta =   0; yo = -1; break;
+            case 1: theta =  45; xo =  1; yo = -1; break;
+            case 2: theta =  90; xo =  1; x = (xsize-1)-x; break;
+            case 3: theta = 135; xo =  1; yo = -1; break;
+            case 4: theta = 180; yo =  1; y = (ysize-1)-y; break;
+            case 5: theta = 225; xo = -1; yo = 1; y = (ysize-1)-y; break;
+            case 6: theta = 270; xo = -1; break;
+            case 7: theta = 315; xo = -1; yo = -1; break;
+            default: return;
+          }
+
+          if (can->getPixelQuick(x,y) != 0
+              || x+xo < 0
+              || x+xo >= xsize
+              || y+yo < 0
+              || y+yo >= ysize )
+          {
+            double gx = feats.get(x, y, z, 1);
+            double gy = feats.get(x, y, z, 2);
+
+            //double go = atan2 (gy, gx);
+
+            norm(x, y)   = sqrt(gx*gx+gy*gy);
+            orient(x, y) = ( gx*cos(theta)+gy*sin(theta) ) / norm(x,y);
+            dist(x, y)   = 0;
+          }
+          else
+          {
+            orient(x, y) = orient(x+xo,y+yo);
+            norm(x, y)   = norm(x+xo,y+yo);
+            dist(x, y)   = dist(x+xo,y+yo) + 1;
+          }
+        }
+
+      for (int y = 0; y < ysize; y++)
+        for (int x = 0; x < xsize; x++)
+        {
+          // distance feature maps
+          feats.set( x, y, z, dist(x,y), firstChannel + dir );
+          // norm feature maps
+          feats.set( x, y, z, norm(x,y), firstChannel + amountDirs + dir );
+          // orientation feature maps
+          feats.set( x, y, z, norm(x,y), firstChannel + (amountDirs*2) + dir );
+        }
+    }
+
+    delete can;
+  }
+}
+
+void SemSegContextTree3D::updateProbabilityMaps (
+    const NICE::MultiChannelImage3DT<unsigned short int> &nodeIndices,
+    NICE::MultiChannelImage3DT<double> &feats,
+    int firstChannel )
+{
+  int xsize = feats.width();
+  int ysize = feats.height();
+  int zsize = feats.depth();
+
+  int classes = ( int ) forest[0][0].dist.size();
+
+  // integral images for context channels (probability maps for each class)
+#pragma omp parallel for
+    for ( int c = 0; c < classes; c++ )
+    {
+      for ( int z = 0; z < zsize; z++ )
+      {
+        for ( int y = 0; y < ysize; y++ )
+        {
+          for ( int x = 0; x < xsize; x++ )
+          {
+            double val = getMeanProb ( x, y, z, c, nodeIndices );
+
+            if (useFeat3)
+              feats ( x, y, z, firstChannel + c ) = val;
+
+            if (useFeat4)
+              feats ( x, y, z, firstChannel + classes + c ) = val;
+          }
+        }
+
+        // Gaussian filter on probability maps
+//        NICE::ImageT<double> img = feats.getChannelT( z, firstChannel+c );
+//        NICE::ImageT<double> gF(xsize,ysize);
+//        NICE::FilterT<double,double,double> filt;
+//        filt.filterGaussSigmaApproximate( img, 2, &gF );
+//        for ( int y = 0; y < ysize; y++ )
+//          for ( int x = 0; x < xsize; x++ )
+//            feats.set(x, y, z, gF.getPixelQuick(x,y), firstChannel+c);
+      }
+
+      feats.calcIntegral ( firstChannel + c );
+    }
+}
+
+inline double computeWeight ( const int &d, const int &dim )
+{
+  if (d == 0)
+    return 0.0;
+  else
+    return 1.0 / ( pow ( 2, ( double ) ( dim - d + 1 ) ) );
+}
+
+void SemSegContextTree3D::train ( const MultiDataset *md )
+{
+  const LabeledSet trainSet = * ( *md ) ["train"];
+  const LabeledSet *trainp = &trainSet;
+  
+  if ( saveLoadData )
+  {
+    if ( FileMgt::fileExists ( fileLocation ) )
+      read ( fileLocation );
+    else
+    {
+      train ( trainp );
+      write ( fileLocation );
+    }
+  }
+  else
+  {
+    train ( trainp );
+  }
+}
+
+void SemSegContextTree3D::train ( const LabeledSet * trainp )
+{
+  int shortsize = numeric_limits<short>::max();
+
+  Timer timer;
+  timer.start();
+  
+  vector<int> zsizeVec;
+  getDepthVector ( trainp, zsizeVec, run3Dseg );
+
+  //FIXME: memory usage
+  vector<MultiChannelImage3DT<double> > allfeats;
+  vector<MultiChannelImage3DT<unsigned short int> > nodeIndices;
+  vector<MultiChannelImageT<int> > labels;
+
+  vector<SparseVector*> globalCategorFeats;
+  vector<map<int,int> > classesPerImage;
+
+  vector<vector<int> > rSize;
+  vector<int> amountRegionpI;
+
+  std::string forbidden_classes_s = conf->gS ( "analysis", "forbidden_classes", "" );
+  classnames.getSelection ( forbidden_classes_s, forbidden_classes );
+
+  int imgCounter = 0;
+  int amountPixels = 0;
+
+  // How many channels of non-integral type do we have?
+
+  if ( imagetype == IMAGETYPE_RGB )
+    rawChannels = 3;
+  else
+    rawChannels = 1;
+
+  if ( useGradient )
+    rawChannels *= 3;
+
+  if ( useWeijer )
+    rawChannels += 11;
+
+  if ( useHoiemFeatures )
+    rawChannels += 8;
+
+  if ( useAdditionalLayer )
+    rawChannels += 1;
+
+
+///////////////////////////// read input data /////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+  int depthCount = 0;
+  vector< string > filelist;
+  NICE::MultiChannelImageT<uchar> pixelLabels;
+
+  for (LabeledSet::const_iterator it = trainp->begin(); it != trainp->end(); it++)
+  {
+    for (std::vector<ImageInfo *>::const_iterator jt = it->second.begin();
+         jt != it->second.end(); jt++)
+    {
+      int classno = it->first;
+      ImageInfo & info = *(*jt);
+      std::string file = info.img();
+      filelist.push_back ( file );
+      depthCount++;
+
+      const LocalizationResult *locResult = info.localization();
+
+      // getting groundtruth
+      NICE::Image pL;
+      pL.resize ( locResult->xsize, locResult->ysize );
+      pL.set ( 0 );
+      locResult->calcLabeledImage ( pL, ( *classNames ).getBackgroundClass() );
+      pixelLabels.addChannel ( pL );
+
+      if ( locResult->size() <= 0 )
+      {
+        fprintf ( stderr, "WARNING: NO ground truth polygons found for %s !\n",
+                  file.c_str() );
+        continue;
+      }
+
+      fprintf ( stderr, "SSContext: Collecting pixel examples from localization info: %s\n", file.c_str() );
+
+      int depthBoundary = 0;
+      if ( run3Dseg )
+      {
+        depthBoundary = zsizeVec[imgCounter];
+      }
+
+      if ( depthCount < depthBoundary ) continue;
+
+      // all image slices collected -> make a 3d image
+      NICE::MultiChannelImage3DT<double> imgData;
+      make3DImage ( filelist, imgData );
+
+      int xsize = imgData.width();
+      int ysize = imgData.height();
+      int zsize = imgData.depth();
+      amountPixels += xsize * ysize * zsize;
+
+      MultiChannelImageT<int> tmpMat ( xsize, ysize, ( uint ) zsize );
+      labels.push_back ( tmpMat );
+
+      nodeIndices.push_back ( MultiChannelImage3DT<unsigned short int> ( xsize, ysize, zsize, nbTrees ) );
+      nodeIndices[imgCounter].setAll ( 0 );
+
+//      MultiChannelImage3DT<double> feats;
+//      allfeats.push_back ( feats );
+
+      int amountRegions;
+      // convert color to L*a*b, add selected feature channels
+      addFeatureMaps ( imgData, filelist, amountRegions );
+      allfeats.push_back(imgData);
+
+      if ( useFeat1 )
+      {
+        amountRegionpI.push_back ( amountRegions );
+        rSize.push_back ( vector<int> ( amountRegions, 0 ) );
+      }
+
+      if ( useCategorization )
+      {
+        globalCategorFeats.push_back ( new SparseVector() );
+        classesPerImage.push_back ( map<int,int>() );
+      }
+
+      for ( int x = 0; x < xsize; x++ )
+      {
+        for ( int y = 0; y < ysize; y++ )
+        {
+          for ( int z = 0; z < zsize; z++ )
+          {
+            if ( useFeat1 )
+              rSize[imgCounter][allfeats[imgCounter] ( x, y, z, rawChannels ) ]++;
+
+            if ( run3Dseg )
+              classno = pixelLabels ( x, y, ( uint ) z );
+            else
+              classno = pL.getPixelQuick ( x,y );
+
+            labels[imgCounter].set ( x, y, classno, ( uint ) z );
+
+            if ( forbidden_classes.find ( classno ) != forbidden_classes.end() )
+              continue;
+
+            labelcounter[classno]++;
+
+            if ( useCategorization )
+              classesPerImage[imgCounter][classno] = 1;
+          }
+        }
+      }
+
+      filelist.clear();
+      pixelLabels.reInit ( 0,0,0 );
+      depthCount = 0;
+      imgCounter++;
+    }
+  }
+
+  int classes = 0;
+  for ( map<int, int>::const_iterator mapit = labelcounter.begin();
+        mapit != labelcounter.end(); mapit++ )
+  {
+    labelmap[mapit->first] = classes;
+    labelmapback[classes] = mapit->first;
+    classes++;
+  }
+
+////////////////////////// channel type configuration /////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+
+  // Type 0: single pixel & pixel-comparison features on gray value channels
+  for ( int i = 0; i < rawChannels; i++ )
+    channelType.push_back ( 0 );
+
+  // Type 1: region channel with unsupervised segmentation
+  int shift = 0;
+  if ( useFeat1 )
+  {
+    channelType.push_back ( 1 );
+    shift++;
+  }
+
+  // Type 2: rectangular and Haar-like features on gray value integral channels
+  if ( useFeat2 )
+    for ( int i = 0; i < rawChannels; i++ )
+      channelType.push_back ( 2 );
+
+  // Type 3: type 2 features on context channels
+  if ( useFeat3 )
+    for ( int i = 0; i < classes; i++ )
+      channelType.push_back ( 3 );
+
+  // Type 4: type 0 features on context channels
+  if ( useFeat4 )
+    for ( int i = 0; i < classes; i++ )
+      channelType.push_back ( 4 );
+
+  // Type 5: ray features for shape modeling on canny-map
+  if ( useFeat5 )
+    for ( int i = 0; i < 24; i++ )
+      channelType.push_back ( 5 );
+
+  // 'amountTypes' sets upper bound for usable feature types
+  int amountTypes = 6;
+  channelsPerType = vector<vector<int> > ( amountTypes, vector<int>() );
+
+  for ( int i = 0; i < ( int ) channelType.size(); i++ )
+  {
+    channelsPerType[channelType[i]].push_back ( i );
+  }
+
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+
+  vector<vector<vector<double> > > regionProbs;
+  if ( useFeat1 )
+  {
+    for ( int i = 0; i < imgCounter; i++ )
+    {
+      regionProbs.push_back ( vector<vector<double> > ( amountRegionpI[i], vector<double> ( classes, 0.0 ) ) );
+    }
+  }
+
+  //balancing
+  a = vector<double> ( classes, 0.0 );
+  int featcounter = 0;
+  for ( int iCounter = 0; iCounter < imgCounter; iCounter++ )
+  {
+    int xsize = ( int ) nodeIndices[iCounter].width();
+    int ysize = ( int ) nodeIndices[iCounter].height();
+    int zsize = ( int ) nodeIndices[iCounter].depth();
+
+    for ( int x = 0; x < xsize; x++ )
+    {
+      for ( int y = 0; y < ysize; y++ )
+      {
+        for ( int z = 0; z < zsize; z++ )
+        {
+          featcounter++;
+          int cn = labels[iCounter] ( x, y, ( uint ) z );
+          if ( labelmap.find ( cn ) == labelmap.end() )
+            continue;
+          a[labelmap[cn]] ++;
+        }
+      }
+    }
+  }
+
+  for ( int i = 0; i < ( int ) a.size(); i++ )
+  {
+    a[i] /= ( double ) featcounter;
+  }
+
+#ifdef VERBOSE
+  cout << "\nDistribution:" << endl;
+  for ( int i = 0; i < ( int ) a.size(); i++ )
+    cout << "class " << i << ": " << a[i] << endl;
+#endif
+
+  depth = 0;
+  uniquenumber = 0;
+
+  //initialize random forest
+  for ( int t = 0; t < nbTrees; t++ )
+  {
+    vector<TreeNode> singletree;
+    singletree.push_back ( TreeNode() );
+    singletree[0].dist = vector<double> ( classes, 0.0 );
+    singletree[0].depth = depth;
+    singletree[0].featcounter = amountPixels;
+    singletree[0].nodeNumber = uniquenumber;
+    uniquenumber++;
+    forest.push_back ( singletree );
+  }
+
+  vector<int> startnode ( nbTrees, 0 );
+
+  bool noNewSplit = false;
+
+  timer.stop();
+  cout << "\nTime for Pre-Processing: " << timer.getLastAbsolute() << " seconds\n" << endl;
+
+  //////////////////////////// train the classifier ///////////////////////////
+  /////////////////////////////////////////////////////////////////////////////
+  timer.start();
+  while ( !noNewSplit && (depth < maxDepth) )
+  {
+    depth++;
+#ifdef DEBUG
+    cout << "depth: " << depth << endl;
+#endif
+    noNewSplit = true;
+    vector<MultiChannelImage3DT<unsigned short int> > lastNodeIndices = nodeIndices;
+    vector<vector<vector<double> > > lastRegionProbs = regionProbs;
+
+    if ( useFeat1 )
+      for ( int i = 0; i < imgCounter; i++ )
+      {
+        int numRegions = (int) regionProbs[i].size();
+        for ( int r = 0; r < numRegions; r++ )
+          for ( int c = 0; c < classes; c++ )
+            regionProbs[i][r][c] = 0.0;
+      }
+
+    // initialize & update context channels
+    for ( int i = 0; i < imgCounter; i++)
+      if ( useFeat3 || useFeat4 )
+        this->updateProbabilityMaps ( nodeIndices[i], allfeats[i], rawChannels + shift );
+
+#ifdef VERBOSE
+    Timer timerDepth;
+    timerDepth.start();
+#endif
+
+    double weight = computeWeight ( depth, maxDepth )
+                    - computeWeight ( depth - 1, maxDepth );
+
+#pragma omp parallel for
+    // for each tree
+    for ( int tree = 0; tree < nbTrees; tree++ )
+    {
+      const int t = ( int ) forest[tree].size();
+      const int s = startnode[tree];
+      startnode[tree] = t;
+      double bestig;
+
+      // for each node
+      for ( int node = s; node < t; node++ )
+      {
+        if ( !forest[tree][node].isleaf && forest[tree][node].left < 0 )
+        {
+          // find best split
+          Operation3D *splitfeat = NULL;
+          double splitval;
+          bestig = getBestSplit ( allfeats, lastNodeIndices, labels, node,
+                                   splitfeat, splitval, tree, lastRegionProbs );
+
+          forest[tree][node].feat = splitfeat;
+          forest[tree][node].decision = splitval;
+
+          // split the node
+          if ( splitfeat != NULL )
+          {
+            noNewSplit = false;
+            int left;
+#pragma omp critical
+            {
+              left = forest[tree].size();
+              forest[tree].push_back ( TreeNode() );
+              forest[tree].push_back ( TreeNode() );
+            }
+            int right = left + 1;
+            forest[tree][node].left = left;
+            forest[tree][node].right = right;
+            forest[tree][left].init( depth, classes, uniquenumber);
+            int leftu = uniquenumber;
+            uniquenumber++;
+            forest[tree][right].init( depth, classes, uniquenumber);
+            int rightu = uniquenumber;
+            uniquenumber++;
+
+#pragma omp parallel for
+            for ( int i = 0; i < imgCounter; i++ )
+            {
+              int xsize = nodeIndices[i].width();
+              int ysize = nodeIndices[i].height();
+              int zsize = nodeIndices[i].depth();
+
+              for ( int x = 0; x < xsize; x++ )
+              {
+                for ( int y = 0; y < ysize; y++ )
+                {
+                  for ( int z = 0; z < zsize; z++ )
+                  {
+                    if ( nodeIndices[i].get ( x, y, z, tree ) == node )
+                    {
+                      // get feature value
+                      Features feat;
+                      feat.feats = &allfeats[i];
+                      feat.rProbs = &lastRegionProbs[i];
+                      double val = 0.0;
+                      val = splitfeat->getVal ( feat, x, y, z );
+                      if ( !isfinite ( val ) ) val = 0.0;
+
+#pragma omp critical
+                      {
+                        int curLabel = labels[i] ( x, y, ( uint ) z );
+                        // traverse to left child
+                        if ( val < splitval )
+                        {
+                          nodeIndices[i].set ( x, y, z, left, tree );
+                          if ( labelmap.find ( curLabel ) != labelmap.end() )
+                            forest[tree][left].dist[labelmap[curLabel]]++;
+                          forest[tree][left].featcounter++;
+                          if ( useCategorization && leftu < shortsize )
+                            ( *globalCategorFeats[i] ) [leftu]+=weight;
+                        }
+                        // traverse to right child
+                        else
+                        {
+                          nodeIndices[i].set ( x, y, z, right, tree );
+                          if ( labelmap.find ( curLabel ) != labelmap.end() )
+                            forest[tree][right].dist[labelmap[curLabel]]++;
+                          forest[tree][right].featcounter++;
+
+                          if ( useCategorization && rightu < shortsize )
+                            ( *globalCategorFeats[i] ) [rightu]+=weight;
+                        }
+                      }
+                    }
+                  }
+                }
+              }
+            }
+
+            // normalize distributions in child leaves
+            double lcounter = 0.0, rcounter = 0.0;
+            for ( int c = 0; c < (int)forest[tree][left].dist.size(); c++ )
+            {
+              if ( forbidden_classes.find ( labelmapback[c] ) != forbidden_classes.end() )
+              {
+                forest[tree][left].dist[c] = 0;
+                forest[tree][right].dist[c] = 0;
+              }
+              else
+              {
+                forest[tree][left].dist[c] /= a[c];
+                lcounter += forest[tree][left].dist[c];
+                forest[tree][right].dist[c] /= a[c];
+                rcounter += forest[tree][right].dist[c];
+              }
+            }
+
+            assert ( lcounter > 0 && rcounter > 0 );
+
+//            if ( lcounter <= 0 || rcounter <= 0 )
+//            {
+//              cout << "lcounter : " << lcounter << " rcounter: " << rcounter << endl;
+//              cout << "splitval: " << splitval << " splittype: " << splitfeat->writeInfos() << endl;
+//              cout << "bestig: " << bestig << endl;
+
+//              for ( int i = 0; i < imgCounter; i++ )
+//              {
+//                int xsize = nodeIndices[i].width();
+//                int ysize = nodeIndices[i].height();
+//                int zsize = nodeIndices[i].depth();
+//                int counter = 0;
+
+//                for ( int x = 0; x < xsize; x++ )
+//                {
+//                  for ( int y = 0; y < ysize; y++ )
+//                  {
+//                    for ( int z = 0; z < zsize; z++ )
+//                    {
+//                      if ( lastNodeIndices[i].get ( x, y, tree ) == node )
+//                      {
+//                        if ( ++counter > 30 )
+//                          break;
+
+//                        Features feat;
+//                        feat.feats = &allfeats[i];
+//                        feat.rProbs = &lastRegionProbs[i];
+
+//                        double val = splitfeat->getVal ( feat, x, y, z );
+//                        if ( !isfinite ( val ) ) val = 0.0;
+
+//                        cout << "splitval: " << splitval << " val: " << val << endl;
+//                      }
+//                    }
+//                  }
+//                }
+//              }
+
+//              assert ( lcounter > 0 && rcounter > 0 );
+//            }
+
+            for ( int c = 0; c < classes; c++ )
+            {
+              forest[tree][left].dist[c] /= lcounter;
+              forest[tree][right].dist[c] /= rcounter;
+            }
+          }
+          else
+          {
+            forest[tree][node].isleaf = true;
+          }
+        }
+      }
+    }
+
+
+    if ( useFeat1 )
+    {
+      for ( int i = 0; i < imgCounter; i++ )
+      {
+        int xsize = nodeIndices[i].width();
+        int ysize = nodeIndices[i].height();
+        int zsize = nodeIndices[i].depth();
+
+#pragma omp parallel for
+        // set region probability distribution
+        for ( int x = 0; x < xsize; x++ )
+        {
+          for ( int y = 0; y < ysize; y++ )
+          {
+            for ( int z = 0; z < zsize; z++ )
+            {
+              for ( int tree = 0; tree < nbTrees; tree++ )
+              {
+                int node = nodeIndices[i].get ( x, y, z, tree );
+                for ( int c = 0; c < classes; c++ )
+                {
+                  int r = (int) ( allfeats[i] ( x, y, z, rawChannels ) );
+                  regionProbs[i][r][c] += forest[tree][node].dist[c];
+                }
+              }
+            }
+          }
+        }
+
+        // normalize distribution
+        int numRegions = (int) regionProbs[i].size();
+        for ( int r = 0; r < numRegions; r++ )
+        {
+          for ( int c = 0; c < classes; c++ )
+          {
+            regionProbs[i][r][c] /= ( double ) ( rSize[i][r] );
+          }
+        }
+      }
+    }
+
+    if ( firstiteration ) firstiteration = false;
+
+#ifdef VERBOSE
+    timerDepth.stop();
+    cout << "Depth " << depth << ": " << timerDepth.getLastAbsolute() << " seconds" <<endl;
+#endif
+
+    lastNodeIndices.clear();
+    lastRegionProbs.clear();
+  }
+
+  timer.stop();
+  cout << "Time for Learning: " << timer.getLastAbsolute() << " seconds\n" << endl;
+
+  //////////////////////// classification using HIK ///////////////////////////
+  /////////////////////////////////////////////////////////////////////////////
+
+  if ( useCategorization && fasthik != NULL )
+  {
+    timer.start();
+    uniquenumber = std::min ( shortsize, uniquenumber );
+    for ( uint i = 0; i < globalCategorFeats.size(); i++ )
+    {
+      globalCategorFeats[i]->setDim ( uniquenumber );
+      globalCategorFeats[i]->normalize();
+    }
+    map<int,Vector> ys;
+
+    int cCounter = 0;
+    for ( map<int,int>::const_iterator it = labelmap.begin();
+          it != labelmap.end(); it++, cCounter++ )
+    {
+      ys[cCounter] = Vector ( globalCategorFeats.size() );
+      for ( int i = 0; i < imgCounter; i++ )
+      {
+        if ( classesPerImage[i].find ( it->first ) != classesPerImage[i].end() )
+        {
+          ys[cCounter][i] = 1;
+        }
+        else
+        {
+          ys[cCounter][i] = -1;
+        }
+      }
+    }
+
+    fasthik->train( reinterpret_cast<vector<const NICE::SparseVector *>&>(globalCategorFeats), ys);
+
+    timer.stop();
+    cerr << "Time for Categorization: " << timer.getLastAbsolute() << " seconds\n" << endl;
+  }
+
+#ifdef VERBOSE
+  cout << "\nFEATURE USAGE" << endl;
+  cout << "#############\n" << endl;
+
+  // amount of used features per feature type
+  std::map<int, int> featTypeCounter;
+  for ( int tree = 0; tree < nbTrees; tree++ )
+  {
+    int t = ( int ) forest[tree].size();
+
+    for ( int node = 0; node < t; node++ )
+    {
+      if ( !forest[tree][node].isleaf && forest[tree][node].left != -1 )
+      {
+        featTypeCounter[ forest[tree][node].feat->getFeatType() ] += 1;
+      }
+    }
+  }
+
+  cout << "Types:" << endl;
+  for ( map<int, int>::const_iterator it = featTypeCounter.begin(); it != featTypeCounter.end(); it++ )
+    cout << it->first << ": " << it->second << endl;
+
+  cout << "\nOperations - All:" << endl;
+  // used operations
+  vector<int> opOverview ( NBOPERATIONS, 0 );
+  // relative use of context vs raw features per tree level
+  vector<vector<double> > contextOverview ( maxDepth, vector<double> ( 2, 0.0 ) );
+  for ( int tree = 0; tree < nbTrees; tree++ )
+  {
+    int t = ( int ) forest[tree].size();
+
+    for ( int node = 0; node < t; node++ )
+    {
+#ifdef DEBUG
+      printf ( "tree[%i]: left: %i, right: %i", node, forest[tree][node].left, forest[tree][node].right );
+#endif
+
+      if ( !forest[tree][node].isleaf && forest[tree][node].left != -1 )
+      {
+        cout <<  forest[tree][node].feat->writeInfos() << endl;
+        opOverview[ forest[tree][node].feat->getOps() ]++;
+        contextOverview[forest[tree][node].depth][ ( int ) forest[tree][node].feat->getContext() ]++;
+      }
+#ifdef DEBUG
+      for ( int d = 0; d < ( int ) forest[tree][node].dist.size(); d++ )
+      {
+        cout << " " << forest[tree][node].dist[d];
+      }
+      cout << endl;
+#endif
+    }
+  }
+
+  // amount of used features per operation type
+  cout << "\nOperations - Summary:" << endl;
+  for ( int t = 0; t < ( int ) opOverview.size(); t++ )
+  {
+    cout << "Ops " << t << ": " << opOverview[ t ] << endl;
+  }
+  // ratio of used context features per depth level
+  cout << "\nContext-Ratio:" << endl;
+  for ( int d = 0; d < maxDepth; d++ )
+  {
+    double sum =  contextOverview[d][0] + contextOverview[d][1];
+    if ( sum == 0 )
+      sum = 1;
+
+    contextOverview[d][0] /= sum;
+    contextOverview[d][1] /= sum;
+
+    cout << "Depth [" << d+1 << "] Normal: " << contextOverview[d][0] << " Context: " << contextOverview[d][1] << endl;
+  }
+#endif
+
+}
+
+void SemSegContextTree3D::addFeatureMaps (
+    NICE::MultiChannelImage3DT<double> &imgData,
+    const vector<string> &filelist,
+    int &amountRegions )
+{
+  int xsize = imgData.width();
+  int ysize = imgData.height();
+  int zsize = imgData.depth();
+
+  amountRegions = 0;
+
+  // RGB to Lab
+  if ( imagetype == IMAGETYPE_RGB )
+  {
+    for ( int z = 0; z < zsize; z++ )
+      for ( int y = 0; y < ysize; y++ )
+        for ( int x = 0; x < xsize; x++ )
+        {
+          double R, G, B, X, Y, Z, L, a, b;
+          R = ( double )imgData.get( x, y, z, 0 ) / 255.0;
+          G = ( double )imgData.get( x, y, z, 1 ) / 255.0;
+          B = ( double )imgData.get( x, y, z, 2 ) / 255.0;
+
+          if ( useAltTristimulus )
+          {
+            ColorConversion::ccRGBtoXYZ( R, G, B, &X, &Y, &Z, 4 );
+            ColorConversion::ccXYZtoCIE_Lab( X, Y, Z, &L, &a, &b, 4 );
+          }
+          else
+          {
+            ColorConversion::ccRGBtoXYZ( R, G, B, &X, &Y, &Z, 0 );
+            ColorConversion::ccXYZtoCIE_Lab( X, Y, Z, &L, &a, &b, 0 );
+          }
+
+          imgData.set( x, y, z, L, 0 );
+          imgData.set( x, y, z, a, 1 );
+          imgData.set( x, y, z, b, 2 );
+        }
+  }
+
+  // Gradient layers
+  if ( useGradient )
+  {
+    int currentsize = imgData.channels();
+    imgData.addChannel ( 2*currentsize );
+
+    for ( int z = 0; z < zsize; z++ )
+      for ( int c = 0; c < currentsize; c++ )
+      {
+        ImageT<double> tmp = imgData.getChannelT(z, c);
+        ImageT<double> sobX( xsize, ysize );
+        ImageT<double> sobY( xsize, ysize );
+        NICE::FilterT<double, double, double>::sobelX ( tmp, sobX );
+        NICE::FilterT<double, double, double>::sobelY ( tmp, sobY );
+        for ( int y = 0; y < ysize; y++ )
+          for ( int x = 0; x < xsize; x++ )
+          {
+            imgData.set( x, y, z, sobX.getPixelQuick(x,y), c+currentsize );
+            imgData.set( x, y, z, sobY.getPixelQuick(x,y), c+(currentsize*2) );
+          }
+      }
+  }
+
+  // Weijer color names
+  if ( useWeijer )
+  {
+    if ( imagetype == IMAGETYPE_RGB )
+    {
+      int currentsize = imgData.channels();
+      imgData.addChannel ( 11 );
+      for ( int z = 0; z < zsize; z++ )
+      {
+        NICE::ColorImage img = imgData.getColor ( z );
+        NICE::MultiChannelImageT<double> cfeats;
+        lfcw->getFeats ( img, cfeats );
+        for ( int c = 0; c < cfeats.channels(); c++)
+          for ( int y = 0; y < ysize; y++ )
+            for ( int x = 0; x < xsize; x++ )
+              imgData.set(x, y, z, cfeats.get(x,y,(uint)c), c+currentsize);
+      }
+    }
+    else
+    {
+      cerr << "Can't compute weijer features of a grayscale image." << endl;
+    }
+  }
+
+  // arbitrary additional layer as image
+  if ( useAdditionalLayer )
+  {
+    int currentsize = imgData.channels();
+    imgData.addChannel ( 1 );
+    for ( int z = 0; z < zsize; z++ )
+    {
+      vector<string> list;
+      StringTools::split ( filelist[z], '/', list );
+      string layerPath = StringTools::trim ( filelist[z], list.back() ) + "addlayer/" + list.back();
+      NICE::Image layer ( layerPath );
+      for ( int y = 0; y < ysize; y++ )
+        for ( int x = 0; x < xsize; x++ )
+          imgData.set(x, y, z, layer.getPixelQuick(x,y), currentsize);
+    }
+  }
+    
+  // read the geometric cues produced by Hoiem et al.
+  if ( useHoiemFeatures )
+  {
+    string hoiemDirectory = conf->gS ( "Features", "hoiem_directory" );
+    // we could also give the following set as a config option
+    string hoiemClasses_s = "sky 000 090-045 090-090 090-135 090 090-por 090-sol";
+    vector<string> hoiemClasses;
+    StringTools::split ( hoiemClasses_s, ' ', hoiemClasses );
+
+    int currentsize = imgData.channels();
+    imgData.addChannel ( hoiemClasses.size() );
+    for ( int z = 0; z < zsize; z++ )
+    {
+      FileName fn ( filelist[z] );
+      fn.removeExtension();
+      FileName fnBase = fn.extractFileName();
+
+      for ( vector<string>::const_iterator i = hoiemClasses.begin(); i != hoiemClasses.end(); i++, currentsize++ )
+      {
+        string hoiemClass = *i;
+        FileName fnConfidenceImage ( hoiemDirectory + fnBase.str() + "_c_" + hoiemClass + ".png" );
+        if ( ! fnConfidenceImage.fileExists() )
+        {
+          fthrow ( Exception, "Unable to read the Hoiem geometric confidence image: " << fnConfidenceImage.str() << " (original image is " << filelist[z] << ")" );
+        }
+        else
+        {
+          Image confidenceImage ( fnConfidenceImage.str() );
+          if ( confidenceImage.width() != xsize || confidenceImage.height() != ysize )
+          {
+            fthrow ( Exception, "The size of the geometric confidence image does not match with the original image size: " << fnConfidenceImage.str() );
+          }
+
+          // copy standard image to double image
+          for ( int y = 0 ; y < confidenceImage.height(); y++ )
+            for ( int x = 0 ; x < confidenceImage.width(); x++ )
+              imgData ( x, y, z, currentsize ) = ( double ) confidenceImage ( x, y );
+
+          currentsize++;
+        }
+      }
+    }
+  }
+
+  // region feature (unsupervised segmentation)
+  int shift = 0;
+  if ( useFeat1 )
+  {
+    shift = 1;
+    MultiChannelImageT<int> regions;
+    regions.reInit( xsize, ysize, zsize );
+    amountRegions = segmentation->segRegions ( imgData, regions, imagetype );
+
+    int currentsize = imgData.channels();
+    imgData.addChannel ( 1 );
+
+    for ( int z = 0; z < ( int ) regions.channels(); z++ )
+      for ( int y = 0; y < regions.height(); y++ )
+        for ( int x = 0; x < regions.width(); x++ )
+          imgData.set ( x, y, z, regions ( x, y, ( uint ) z ), currentsize );
+
+  }
+
+  // intergal images of raw channels
+  if ( useFeat2 )
+  {
+    imgData.addChannel ( rawChannels );
+
+#pragma omp parallel for
+    for ( int i = 0; i < rawChannels; i++ )
+    {
+      int corg = i;
+      int cint = i + rawChannels + shift;
+
+      for ( int z = 0; z < zsize; z++ )
+        for ( int y = 0; y < ysize; y++ )
+          for ( int x = 0; x < xsize; x++ )
+            imgData ( x, y, z, cint ) = imgData ( x, y, z, corg );
+
+      imgData.calcIntegral ( cint );
+    }
+  }
+
+  int classes = classNames->numClasses();
+
+  if ( useFeat3 )
+    imgData.addChannel ( classes );
+
+  if ( useFeat4 )
+    imgData.addChannel ( classes );
+
+  if ( useFeat5 )
+  {
+    imgData.addChannel ( 24 );
+    this->computeRayFeatImage( imgData, imgData.channels()-24);
+  }
+
+}
+
+void SemSegContextTree3D::classify (
+    NICE::MultiChannelImage3DT<double> & imgData,
+    NICE::MultiChannelImageT<double> & segresult,
+    NICE::MultiChannelImage3DT<double> & probabilities,
+    const std::vector<std::string> & filelist )
+{
+  int xsize = imgData.width();
+  int ysize = imgData.height();
+  int zsize = imgData.depth();
+
+  ////////////////////////// initialize variables /////////////////////////////
+  /////////////////////////////////////////////////////////////////////////////
+
+  firstiteration = true;
+  depth = 0;
+
+  Timer timer;
+  timer.start();
+
+  // classes occurred during training step
+  int classes = labelmapback.size();
+  // classes defined in config file
+  int numClasses = classNames->numClasses();
+
+  // class probabilities by pixel
+  probabilities.reInit ( xsize, ysize, zsize, numClasses );
+  probabilities.setAll ( 0 );
+
+  // class probabilities by region
+  vector<vector<double> > regionProbs;
+
+  // affiliation: pixel <-> (tree,node)
+  MultiChannelImage3DT<unsigned short int> nodeIndices ( xsize, ysize, zsize, nbTrees );
+  nodeIndices.setAll ( 0 );
+
+  // for categorization
+  SparseVector *globalCategorFeat;
+  globalCategorFeat = new SparseVector();
+
+  /////////////////////////// get feature values //////////////////////////////
+  /////////////////////////////////////////////////////////////////////////////
+
+  // Basic Features
+  int amountRegions;
+  addFeatureMaps ( imgData, filelist, amountRegions );
+
+  vector<int> rSize;
+  int shift = 0;
+  if ( useFeat1 )
+  {
+    shift = 1;
+    regionProbs = vector<vector<double> > ( amountRegions, vector<double> ( classes, 0.0 ) );
+    rSize = vector<int> ( amountRegions, 0 );
+    for ( int z = 0; z < zsize; z++ )
+    {
+      for ( int y = 0; y < ysize; y++ )
+      {
+        for ( int x = 0; x < xsize; x++ )
+        {
+          rSize[imgData ( x, y, z, rawChannels ) ]++;
+        }
+      }
+    }
+  }
+
+  ////////////////// traverse image example through trees /////////////////////
+  /////////////////////////////////////////////////////////////////////////////
+
+  bool noNewSplit = false;
+  for ( int d = 0; d < maxDepth && !noNewSplit; d++ )
+  {
+    depth++;
+    vector<vector<double> > lastRegionProbs = regionProbs;
+
+    if ( useFeat1 )
+    {
+      int numRegions = ( int ) regionProbs.size();
+      for ( int r = 0; r < numRegions; r++ )
+        for ( int c = 0; c < classes; c++ )
+          regionProbs[r][c] = 0.0;
+    }
+
+    if ( depth < maxDepth )
+    {
+      int firstChannel = rawChannels + shift;
+      if ( useFeat3 || useFeat4 )
+        this->updateProbabilityMaps ( nodeIndices, imgData, firstChannel );
+    }
+
+    double weight = computeWeight ( depth, maxDepth )
+                    - computeWeight ( depth - 1, maxDepth );
+
+    noNewSplit = true;
+
+    int tree;
+#pragma omp parallel for private(tree)
+    for ( tree = 0; tree < nbTrees; tree++ )
+    {
+      for ( int x = 0; x < xsize; x++ )
+      {
+        for ( int y = 0; y < ysize; y++ )
+        {
+          for ( int z = 0; z < zsize; z++ )
+          {
+            int node = nodeIndices.get ( x, y, z, tree );
+
+            if ( forest[tree][node].left > 0 )
+            {
+              noNewSplit = false;
+              Features feat;
+              feat.feats = &imgData;
+              feat.rProbs = &lastRegionProbs;
+
+              double val = forest[tree][node].feat->getVal ( feat, x, y, z );
+              if ( !isfinite ( val ) ) val = 0.0;
+
+              // traverse to left child
+              if ( val < forest[tree][node].decision )
+              {
+                int left = forest[tree][node].left;
+                nodeIndices.set ( x, y, z, left, tree );
+#pragma omp critical
+                {
+                  if ( fasthik != NULL
+                       && useCategorization
+                       && forest[tree][left].nodeNumber < uniquenumber )
+                    ( *globalCategorFeat ) [forest[tree][left].nodeNumber] += weight;
+                }
+              }
+              // traverse to right child
+              else
+              {
+                int right = forest[tree][node].right;
+                nodeIndices.set ( x, y, z, right, tree );
+#pragma omp critical
+                {
+                  if ( fasthik != NULL
+                       && useCategorization
+                       && forest[tree][right].nodeNumber < uniquenumber )
+                    ( *globalCategorFeat ) [forest[tree][right].nodeNumber] += weight;
+                }
+              }
+            }
+          }
+        }
+      }
+    }
+
+    if ( useFeat1 )
+    {
+      int xsize = nodeIndices.width();
+      int ysize = nodeIndices.height();
+      int zsize = nodeIndices.depth();
+
+#pragma omp parallel for
+      for ( int x = 0; x < xsize; x++ )
+      {
+        for ( int y = 0; y < ysize; y++ )
+        {
+          for ( int z = 0; z < zsize; z++ )
+          {
+            for ( int tree = 0; tree < nbTrees; tree++ )
+            {
+              int node = nodeIndices.get ( x, y, z, tree );
+              for ( uint c = 0; c < forest[tree][node].dist.size(); c++ )
+              {
+                int r = (int) imgData ( x, y, z, rawChannels );
+                regionProbs[r][c] += forest[tree][node].dist[c];
+              }
+            }
+          }
+        }
+      }
+
+      int numRegions = (int) regionProbs.size();
+      for ( int r = 0; r < numRegions; r++ )
+      {
+        for ( int c = 0; c < (int) classes; c++ )
+        {
+          regionProbs[r][c] /= ( double ) ( rSize[r] );
+        }
+      }
+    }
+
+    if ( (depth < maxDepth) && firstiteration ) firstiteration = false;
+  }
+
+  vector<int> classesInImg;
+
+  if ( useCategorization )
+  {
+    if ( cndir != "" )
+    {
+      for ( int z = 0; z < zsize; z++ )
+      {
+        vector< string > list;
+        StringTools::split ( filelist[z], '/', list );
+        string orgname = list.back();
+
+        ifstream infile ( ( cndir + "/" + orgname + ".dat" ).c_str() );
+        while ( !infile.eof() && infile.good() )
+        {
+          int tmp;
+          infile >> tmp;
+          assert ( tmp >= 0 && tmp < numClasses );
+          classesInImg.push_back ( tmp );
+        }
+      }
+    }
+    else
+    {
+      globalCategorFeat->setDim ( uniquenumber );
+      globalCategorFeat->normalize();
+      ClassificationResult cr = fasthik->classify( globalCategorFeat);
+      for ( uint i = 0; i < ( uint ) classes; i++ )
+      {
+        cerr << cr.scores[i] << " ";
+        if ( cr.scores[i] > 0.0/*-0.3*/ )
+        {
+          classesInImg.push_back ( i );
+        }
+      }
+    }
+    cerr << "amount of classes: " << classes << " used classes: " << classesInImg.size() << endl;
+  }
+
+  if ( classesInImg.size() == 0 )
+  {
+    for ( uint i = 0; i < ( uint ) classes; i++ )
+    {
+      classesInImg.push_back ( i );
+    }
+  }
+
+  // final labeling step
+  if ( pixelWiseLabeling )
+  {
+    for ( int x = 0; x < xsize; x++ )
+    {
+      for ( int y = 0; y < ysize; y++ )
+      {
+        for ( int z = 0; z < zsize; z++ )
+        {
+          //TODO by nodes instead of pixel?
+          double maxProb = - numeric_limits<double>::max();
+          int maxClass = 0;
+
+          for ( uint c = 0; c < classesInImg.size(); c++ )
+          {
+            int i = classesInImg[c];
+            double curProb = getMeanProb ( x, y, z, i, nodeIndices );
+            probabilities ( x, y, z, labelmapback[i] ) = curProb;
+
+            if ( curProb > maxProb )
+            {
+              maxProb = curProb;
+              maxClass = labelmapback[i];
+            }
+          }
+          assert(maxProb <= 1);
+
+          // copy pixel labeling into segresults (output)
+          segresult.set ( x, y, maxClass, ( uint ) z );
+        }
+      }
+    }
+
+#ifdef VISUALIZE
+    getProbabilityMap( probabilities );
+#endif
+  }
+  else
+  {
+    // labeling by region
+    NICE::MultiChannelImageT<int> regions;
+    int xsize = imgData.width();
+    int ysize = imgData.height();
+    int zsize = imgData.depth();
+    regions.reInit ( xsize, ysize, zsize );
+
+    if ( useFeat1 )
+    {
+      int rchannel = -1;
+      for ( uint i = 0; i < channelType.size(); i++ )
+      {
+        if ( channelType[i] == 1 )
+        {
+          rchannel = i;
+          break;
+        }
+      }
+
+      assert ( rchannel > -1 );
+
+      for ( int z = 0; z < zsize; z++ )
+      {
+        for ( int y = 0; y < ysize; y++ )
+        {
+          for ( int x = 0; x < xsize; x++ )
+          {
+            regions.set ( x, y, imgData ( x, y, z, rchannel ), ( uint ) z );
+          }
+        }
+      }
+    }
+    else
+    {
+      amountRegions = segmentation->segRegions ( imgData, regions, imagetype );
+
+#ifdef DEBUG
+      for ( unsigned int z = 0; z < ( uint ) zsize; z++ )
+      {
+        NICE::Matrix regmask;
+        NICE::ColorImage colorimg ( xsize, ysize );
+        NICE::ColorImage marked ( xsize, ysize );
+        regmask.resize ( xsize, ysize );
+        for ( int y = 0; y < ysize; y++ )
+        {
+          for ( int x = 0; x < xsize; x++ )
+          {
+            regmask ( x,y ) = regions ( x,y,z );
+            colorimg.setPixelQuick ( x, y, 0, imgData.get ( x,y,z,0 ) );
+            colorimg.setPixelQuick ( x, y, 1, imgData.get ( x,y,z,0 ) );
+            colorimg.setPixelQuick ( x, y, 2, imgData.get ( x,y,z,0 ) );
+          }
+        }
+        vector<int> colorvals;
+        colorvals.push_back ( 255 );
+        colorvals.push_back ( 0 );
+        colorvals.push_back ( 0 );
+        segmentation->markContours ( colorimg, regmask, colorvals, marked );
+        std::vector<string> list;
+        StringTools::split ( filelist[z], '/', list );
+        string savePath = StringTools::trim ( filelist[z], list.back() ) + "marked/" + list.back();
+        marked.write ( savePath );
+      }
+#endif
+    }
+
+    regionProbs.clear();
+    regionProbs = vector<vector<double> > ( amountRegions, vector<double> ( classes, 0.0 ) );
+
+    vector<int> bestlabels ( amountRegions, labelmapback[classesInImg[0]] );
+    for ( int z = 0; z < zsize; z++ )
+    {
+      for ( int y = 0; y < ysize; y++ )
+      {
+        for ( int x = 0; x < xsize; x++ )
+        {
+          int r = regions ( x, y, ( uint ) z );
+          for ( uint i = 0; i < classesInImg.size(); i++ )
+          {
+            int c = classesInImg[i];
+            // get mean voting of all trees
+            regionProbs[r][c] += getMeanProb ( x, y, z, c, nodeIndices );
+          }
+        }
+      }
+    }
+
+    for ( int r = 0; r < amountRegions; r++ )
+    {
+      double maxProb = regionProbs[r][classesInImg[0]];
+      bestlabels[r] = classesInImg[0];
+
+      for ( int c = 1; c < classes; c++ )
+      {
+        if ( maxProb < regionProbs[r][c] )
+        {
+          maxProb = regionProbs[r][c];
+          bestlabels[r] = c;
+        }
+      }
+
+      bestlabels[r] = labelmapback[bestlabels[r]];
+    }
+
+    // copy region labeling into segresults (output)
+    for ( int z = 0; z < zsize; z++ )
+    {
+      for ( int y = 0; y < ysize; y++ )
+      {
+        for ( int x = 0; x < xsize; x++ )
+        {
+          segresult.set ( x, y, bestlabels[regions ( x,y, ( uint ) z ) ], ( uint ) z );
+        }
+      }
+    }
+
+#ifdef WRITEREGIONS
+    for ( int z = 0; z < zsize; z++ )
+    {
+      RegionGraph rg;
+      NICE::ColorImage img ( xsize,ysize );
+      if ( imagetype == IMAGETYPE_RGB )
+      {
+        img = imgData.getColor ( z );
+      }
+      else
+      {
+        NICE::Image gray = imgData.getChannel ( z );
+        for ( int y = 0; y < ysize; y++ )
+        {
+          for ( int x = 0; x < xsize; x++ )
+          {
+            int val = gray.getPixelQuick ( x,y );
+            img.setPixelQuick ( x, y, val, val, val );
+          }
+        }
+      }
+
+      Matrix regions_tmp ( xsize,ysize );
+      for ( int y = 0; y < ysize; y++ )
+      {
+        for ( int x = 0; x < xsize; x++ )
+        {
+          regions_tmp ( x,y ) = regions ( x,y, ( uint ) z );
+        }
+      }
+      segmentation->getGraphRepresentation ( img, regions_tmp,  rg );
+      for ( uint pos = 0; pos < regionProbs.size(); pos++ )
+      {
+        rg[pos]->setProbs ( regionProbs[pos] );
+      }
+
+      std::string s;
+      std::stringstream out;
+      std::vector< std::string > list;
+      StringTools::split ( filelist[z], '/', list );
+
+      out << "rgout/" << list.back() << ".graph";
+      string writefile = out.str();
+      rg.write ( writefile );
+    }
+#endif
+  }
+
+  timer.stop();
+  cout << "\nTime for Classification: " << timer.getLastAbsolute() << endl;
+
+  // CLEANING UP
+  // TODO: operations in "forest"
+  while( !ops.empty() )
+  {
+    vector<Operation3D*> &tops = ops.back();
+    while ( !tops.empty() )
+      tops.pop_back();
+
+    ops.pop_back();
+  }
+
+  delete globalCategorFeat;
+}
+
+void SemSegContextTree3D::store ( std::ostream & os, int format ) const
+{
+  os.precision ( numeric_limits<double>::digits10 + 1 );
+  os << nbTrees << endl;
+  classnames.store ( os );
+
+  map<int, int>::const_iterator it;
+
+  os << labelmap.size() << endl;
+  for ( it = labelmap.begin() ; it != labelmap.end(); it++ )
+    os << ( *it ).first << " " << ( *it ).second << endl;
+
+  os << labelmapback.size() << endl;
+  for ( it = labelmapback.begin() ; it != labelmapback.end(); it++ )
+    os << ( *it ).first << " " << ( *it ).second << endl;
+
+  int trees = forest.size();
+  os << trees << endl;
+
+  for ( int t = 0; t < trees; t++ )
+  {
+    int nodes = forest[t].size();
+    os << nodes << endl;
+    for ( int n = 0; n < nodes; n++ )
+    {
+      os << forest[t][n].left << " " << forest[t][n].right << " " << forest[t][n].decision << " " << forest[t][n].isleaf << " " << forest[t][n].depth << " " << forest[t][n].featcounter << " " << forest[t][n].nodeNumber << endl;
+      os << forest[t][n].dist << endl;
+
+      if ( forest[t][n].feat == NULL )
+        os << -1 << endl;
+      else
+      {
+        os << forest[t][n].feat->getOps() << endl;
+        forest[t][n].feat->store ( os );
+      }
+    }
+  }
+
+  os << channelType.size() << endl;
+  for ( int i = 0; i < ( int ) channelType.size(); i++ )
+  {
+    os << channelType[i] << " ";
+  }
+  os << endl;
+
+  os << rawChannels << endl;
+
+  os << uniquenumber << endl;
+}
+
+void SemSegContextTree3D::restore ( std::istream & is, int format )
+{
+  is >> nbTrees;
+
+  classnames.restore ( is );
+
+  int lsize;
+  is >> lsize;
+
+  labelmap.clear();
+  for ( int l = 0; l < lsize; l++ )
+  {
+    int first, second;
+    is >> first;
+    is >> second;
+    labelmap[first] = second;
+  }
+
+  is >> lsize;
+  labelmapback.clear();
+  for ( int l = 0; l < lsize; l++ )
+  {
+    int first, second;
+    is >> first;
+    is >> second;
+    labelmapback[first] = second;
+  }
+
+  int trees;
+  is >> trees;
+  forest.clear();
+
+  for ( int t = 0; t < trees; t++ )
+  {
+    vector<TreeNode> tmptree;
+    forest.push_back ( tmptree );
+    int nodes;
+    is >> nodes;
+
+    for ( int n = 0; n < nodes; n++ )
+    {
+      TreeNode tmpnode;
+      forest[t].push_back ( tmpnode );
+      is >> forest[t][n].left;
+      is >> forest[t][n].right;
+      is >> forest[t][n].decision;
+      is >> forest[t][n].isleaf;
+      is >> forest[t][n].depth;
+      is >> forest[t][n].featcounter;
+      is >> forest[t][n].nodeNumber;
+
+      is >> forest[t][n].dist;
+
+      int feattype;
+      is >> feattype;
+      assert ( feattype < NBOPERATIONS );
+      forest[t][n].feat = NULL;
+
+      if ( feattype >= 0 )
+      {
+        for ( uint o = 0; o < ops.size(); o++ )
+        {
+          for ( uint o2 = 0; o2 < ops[o].size(); o2++ )
+          {
+            if ( forest[t][n].feat == NULL )
+            {
+              if ( ops[o][o2]->getOps() == feattype )
+              {
+                forest[t][n].feat = ops[o][o2]->clone();
+                break;
+              }
+            }
+          }
+        }
+
+        assert ( forest[t][n].feat != NULL );
+        forest[t][n].feat->restore ( is );
+
+      }
+    }
+  }
+
+  channelType.clear();
+  int ctsize;
+  is >> ctsize;
+  for ( int i = 0; i < ctsize; i++ )
+  {
+    int tmp;
+    is >> tmp;
+    switch (tmp)
+    {
+      case 0: useFeat0 = true; break;
+      case 1: useFeat1 = true; break;
+      case 2: useFeat2 = true; break;
+      case 3: useFeat3 = true; break;
+      case 4: useFeat4 = true; break;
+      case 5: useFeat5 = true; break;
+    }
+    channelType.push_back ( tmp );
+  }
+
+  // integralMap is deprecated but kept in RESTORE
+  // for downwards compatibility!
+//  std::vector<std::pair<int, int> > integralMap;
+//  integralMap.clear();
+//  int iMapSize;
+//  is >> iMapSize;
+//  for ( int i = 0; i < iMapSize; i++ )
+//  {
+//    int first;
+//    int second;
+//    is >> first;
+//    is >> second;
+//    integralMap.push_back ( pair<int, int> ( first, second ) );
+//  }
+
+  is >> rawChannels;
+
+  is >> uniquenumber;
+}
+
+//###################### DEPRECATED #########################//
+
+void SemSegContextTree3D::semanticseg ( OBJREC::CachedExample *ce,
+                               NICE::Image & segresult,
+                               NICE::MultiChannelImageT<double> & probabilities )
+{}

semseg/SemSegContextTree3D.h

@@ -0,0 +1,285 @@
+/**
+* @file SemSegContextTree3D.h
+* @brief Context Trees -> Combination of decision tree and context information
+* @author Björn Fröhlich, Sven Sickert
+* @date 29.11.2011
+
+*/
+#ifndef SemSegContextTree3DINCLUDE
+#define SemSegContextTree3DINCLUDE
+
+// nice-core includes
+#include <core/vector/VVector.h>
+#include <core/image/MultiChannelImage3DT.h>
+
+// nice-gphik-exp includes
+#include <gp-hik-exp/GPHIKClassifierNICE.h>
+
+// nice-vislearning includes
+#include <vislearning/features/localfeatures/LocalFeatureColorWeijer.h>
+
+// nice-segmentation includes
+#include <segmentation/RegionSegmentationMethod.h>
+
+// nice-semseg includes
+#include "SemanticSegmentation.h"
+#include "semseg3d/semseg/operations/Operations3D.h"
+
+namespace OBJREC
+{
+
+/** Localization system */
+
+class SemSegContextTree3D : public SemanticSegmentation
+{
+private:
+  /** Segmentation Method */
+  RegionSegmentationMethod *segmentation;
+
+  /** tree -> saved as vector of nodes */
+  std::vector<std::vector<TreeNode> > forest;
+
+  /** local features */
+  LocalFeatureColorWeijer *lfcw;
+
+  /** whether to run in 3D mode or not */
+  bool run3Dseg;
+
+  /** whether to use a particular feature type or not */
+  bool useFeat0, useFeat1, useFeat2, useFeat3, useFeat4, useFeat5;
+
+  /** array of usable feature types*/
+  std::vector<int> featTypes;
+
+  /** maximum samples for tree  */
+  int maxSamples;
+
+  /** size for neighbourhood */
+  int windowSize;
+
+  /** multiplier for window size if context feature */
+  int contextMultiplier;
+
+  /** how many feats should be considered for a split */
+  int featsPerSplit;
+
+  /** count samples per label */
+  std::map<int, int> labelcounter;
+
+  /** map of labels */
+  std::map<int, int> labelmap;
+
+  /** map of labels inverse*/
+  std::map<int, int> labelmapback;
+
+  /** scalefactor for balancing for each class */
+  std::vector<double> a;
+
+  /** the minimum number of features allowed in a leaf */
+  int minFeats;
+
+  /** maximal depth of tree */
+  int maxDepth;
+
+  /** current depth for training */
+  int depth;
+
+  /** how many splittests */
+  int randomTests;
+
+  /** prototype operations for pairwise features */
+  std::vector<std::vector<Operation3D*> > ops;
+
+  /** use alternative calculation for information gain */
+  bool useShannonEntropy;
+
+  /** Classnames */
+  ClassNames classnames;
+
+  /** train selection */
+  std::set<int> forbidden_classes;
+
+  /** Configfile */
+  const NICE::Config *conf;
+
+  /** use pixelwise labeling or regionlabeling with additional segmenation */
+  bool pixelWiseLabeling;
+
+  /** Number of trees used for the forest */
+  int nbTrees;
+
+  /** whether to use alternative tristimulus for CIE_Lab that matches openCV or not */
+  bool useAltTristimulus;
+
+  /** use Gradient image or not */
+  bool useGradient;
+
+  /** use Color features from van de Weijer or not */
+  bool useWeijer;
+
+  /** use additional input Layer or not */
+  bool useAdditionalLayer;
+
+  /** use external image categorization to avoid some classes */
+  bool useCategorization;
+
+  /** categorization information for external categorization */
+  std::string cndir;
+
+  /** how to handle each channel
+   * 0: simple grayvalue features
+   * 1: which pixel belongs to which region
+   * 2: grayvalue integral images
+   * 3: context integral images
+   * 4: simple context features
+   */
+  std::vector<int> channelType;
+
+  /** list of channels per feature type */
+  std::vector<std::vector<int> > channelsPerType;
+
+  /** whether we should use the geometric features of Hoeim (only offline computation with MATLAB supported) */
+  bool useHoiemFeatures;
+  
+  /** save / load trained icf classifier */
+  bool saveLoadData;
+
+  /** file location of trained icf classifier */
+  std::string fileLocation;
+
+  /** first iteration or not */
+  bool firstiteration;
+
+  /** amount of grayvalue Channels */
+  int rawChannels;
+
+  /** classifier for categorization */
+  OBJREC::GPHIKClassifierNICE *fasthik;
+
+  /** unique numbers for nodes */
+  int uniquenumber;
+  
+  /**
+   * @brief initOperations initialize the operation types
+   */
+  void initOperations();
+
+  /**
+   * @brief train the actual training method
+   * @param trainp pointer to training data
+   */
+  void train ( const LabeledSet * trainp );
+
+  /**
+   * @brief updateProbabilityMaps computes probability maps for context features
+   * @param nodeIndices matrix with current node for each feature
+   * @param feats output MCI3D (must be initilized)
+   * @param firstChannel index of the first channel
+   */
+  void updateProbabilityMaps ( const NICE::MultiChannelImage3DT<unsigned short int> &nodeIndices, NICE::MultiChannelImage3DT<double> &feats, int firstChannel );
+
+  /**
+   * @brief computeRayFeatImage computes ray feature images using canny filter
+   * @param feats output MCI3D (must be initilized)
+   * @param firstChannel index of the first channel
+   */
+  void computeRayFeatImage ( NICE::MultiChannelImage3DT<double> &feats, int firstChannel );
+
+  /**
+   * @brief addFeatureMaps initializes the selected feature channels
+   * @param imgData output MCI3D (must be initilized)
+   * @param filelist a list of image file names representing slices of a stack
+   * @param amountRegions the amount of regions created by the segmentation
+   **/
+  void addFeatureMaps ( NICE::MultiChannelImage3DT<double> &imgData, const std::vector<std::string> &filelist, int &amountRegions );
+
+  /**
+   * compute best split for current settings
+   * @param feats features
+   * @param nodeIndices matrix with current node for each feature
+   * @param labels labels for each feature
+   * @param node current node
+   * @param splitfeat output selected feature dimension
+   * @param splitval output threshold for selected feature
+   * @return double best information gain value
+   */
+  double getBestSplit ( std::vector<NICE::MultiChannelImage3DT<double> > &feats, std::vector<NICE::MultiChannelImage3DT<unsigned short int> > &nodeIndices, const std::vector<NICE::MultiChannelImageT<int> > &labels, int node, Operation3D *&splitop, double &splitval, const int &tree, std::vector<std::vector<std::vector<double> > > &regionProbs );
+
+  /**
+   * @brief computes the mean probability for a given class over all trees
+   * @param x x position
+   * @param y y position
+   * @param z z position
+   * @param channel current class
+   * @param nodeIndices matrix with current node for each feature
+   * @return double mean value
+   **/
+  inline double getMeanProb ( const int &x, const int &y, const int &z, const int &channel, const NICE::MultiChannelImage3DT<unsigned short int> &nodeIndices );
+
+public:
+  /** simple constructor */
+  SemSegContextTree3D ();
+
+  /** constructor */
+  SemSegContextTree3D ( const NICE::Config *conf, const MultiDataset *md );
+
+  /** simple destructor */
+  virtual ~SemSegContextTree3D();
+
+  /**
+   * classify each pixel of a single 3d image
+   * @param imgData input data
+   * @param segresult segmentation results
+   * @param probabilities probabilities for each pixel
+   */
+  void classify ( NICE::MultiChannelImage3DT<double> &imgData,
+                  NICE::MultiChannelImageT<double> & segresult,
+                  NICE::MultiChannelImage3DT<double> & probabilities,
+                  const std::vector<std::string> & filelist );
+
+  /**
+   * the training method with checking for already existing trained classifier from file
+   * @param md training data
+   */
+  void train ( const MultiDataset *md );
+
+  // deprecated stuff
+  virtual void semanticseg ( OBJREC::CachedExample *ce,
+                             NICE::Image & segresult,
+                             NICE::MultiChannelImageT<double> & probabilities );
+
+  bool active3DMode ()
+  {
+    return run3Dseg;
+  }
+
+  /**
+   * @brief load all data to is stream
+   *
+   * @param is input stream
+   * @param format has no influence
+   * @return void
+   **/
+  virtual void restore ( std::istream & is, int format = 0 );
+
+  /**
+   * @brief save all data to is stream
+   *
+   * @param os output stream
+   * @param format has no influence
+   * @return void
+   **/
+  virtual void store ( std::ostream & os, int format = 0 ) const;
+
+  /**
+   * @brief clean up
+   *
+   * @return void
+   **/
+  virtual void clear () {}
+
+};
+
+} // namespace
+
+#endif

semseg/SemanticSegmentation.cpp

@@ -1,12 +1,15 @@
 /**
 * @file SemanticSegmentation.cpp
 * @brief abstract interface for semantic segmentation algorithms
-* @author Erik Rodner, Alexander Freytag
-* @date 03/19/2009, latest update: 29-01-2014 (dd-mm-yyyy)
+* @author Erik Rodner, Alexander Freytag, Sven Sickert
+* @date 03/19/2009
 
 */
 #include <iostream>
 
+#include "core/image/MultiChannelImage3DT.h"
+#include "core/basics/StringTools.h"
+
 #include "SemanticSegmentation.h"
 #include "vislearning/baselib/Preprocess.h"
 #include "vislearning/baselib/Globals.h"
@@ -15,7 +18,6 @@ using namespace OBJREC;
 using namespace std;
 using namespace NICE;
 
-
 ///////////////////// ///////////////////// /////////////////////
 //                   CONSTRUCTORS / DESTRUCTORS
 ///////////////////// ///////////////////// /////////////////////     
@@ -29,7 +31,6 @@ SemanticSegmentation::SemanticSegmentation ( )
 
 }
 
-
 SemanticSegmentation::SemanticSegmentation ( const Config *conf,
                                              const ClassNames *classNames )
     : iterationCountSuffix(1)
@@ -70,7 +71,7 @@ void SemanticSegmentation::initFromConfig(const Config* conf, const string& s_co
 
 ///////////////////// ///////////////////// /////////////////////
 //                      SEGMENTATION STUFF
-///////////////////// ///////////////////// ///////////////////// 
+///////////////////// ///////////////////// /////////////////////
 
 void SemanticSegmentation::semanticseg ( const std::string & filename,
     NICE::Image & segresult,
@@ -95,7 +96,6 @@ void SemanticSegmentation::semanticseg ( const std::string & filename,
 ///////////////////// ///////////////////// /////////////////////
 //                      DATA CONVERSION
 ///////////////////// ///////////////////// ///////////////////// 
-
 void SemanticSegmentation::convertLSetToSparseExamples ( Examples &examples, LabeledSetVector &lvec )
 {
 #ifdef DEBUG_PRINTS
@@ -140,7 +140,7 @@ void SemanticSegmentation::convertLSetToExamples ( Examples &examples, LabeledSe
       examples.push_back ( pair<int, Example> ( iter->first, ex ) );
     }
   }
-  
+
   if (!removeOldDataPointer)
   {
     //NOTE this is only useful, if our classifier does NOT need the data explicitely
@@ -236,7 +236,6 @@ void SemanticSegmentation::convertVVectorToExamples ( VVector &feats, Examples &
 #endif
 }
 
-
 void SemanticSegmentation::setIterationCountSuffix( const int & _iterationCountSuffix)
 {
   this->iterationCountSuffix = _iterationCountSuffix;
@@ -247,77 +246,224 @@ void SemanticSegmentation::setClassNames ( const OBJREC::ClassNames * _className
   this->classNames = _classNames;
 }
 
+void SemanticSegmentation::getDepthVector ( const LabeledSet *Files,
+                                            vector<int> & depthVec,
+                                            const bool run3Dseg )
+{
+  std::string oldName;
+  int zsize = 0;
+  bool isInit = false;
+
+  for (LabeledSet::const_iterator it = Files->begin(); it != Files->end(); it++)
+  {
+    for (std::vector<ImageInfo *>::const_iterator jt = it->second.begin();
+         jt != it->second.end(); jt++)
+    {
+      ImageInfo & info = *(*jt);
+      std::string file = info.img();
+
+      std::vector< std::string > list;
+      StringTools::split ( file, '/', list );
+      std::string filename = list.back();
+      uint found = filename.find_last_of ( "_" );
+      if (run3Dseg && found < filename.size() && found-3 > 0 )
+      {
+        std::string curName = filename.substr ( found-3,3 );
+        if ( !isInit )
+        {
+          oldName = curName;
+          isInit = true;
+        }
+        if ( curName.compare ( oldName ) == 0 ) // if strings match up
+        {
+          zsize++;
+        }
+        else
+        {
+          depthVec.push_back ( zsize );
+          zsize = 1;
+          oldName = curName;
+        }
+      }
+      else
+      {
+        zsize = 1;
+        depthVec.push_back ( zsize );
+      }
+
+    }
+  }
+  depthVec.push_back ( zsize );
+}
+
+void SemanticSegmentation::make3DImage ( const std::vector<std::string> & filelist,
+    NICE::MultiChannelImage3DT<double> & imgData )
+{
+  bool isInit = false;
+  for ( int it = 0; it < ( int ) filelist.size(); it++ )
+  {
+    if ( imagetype == IMAGETYPE_RGB )
+    {
+      NICE::ColorImage img;
+      try
+      {
+        img.read ( filelist[it] );
+      }
+      catch ( ImageException iE )
+      {
+        fprintf ( stderr, "Failed to open color image file: %s\n", filelist[it].c_str() );
+        fprintf ( stderr, "%s\n", iE.what() );
+        exit ( -1 );
+      }
+      if ( !isInit )
+      {
+        imgData.reInit ( img.width(),img.height(),filelist.size(),3 );
+        isInit = true;
+      }
+      for ( int y = 0; y < img.height(); y++ )
+      {
+        for ( int x = 0; x < img.width(); x++ )
+        {
+          for ( int r = 0; r < 3; r++ )
+          {
+            imgData.set ( x, y, it, img.getPixel ( x,y,r ), r );
+          }
+        }
+      }
+    }
+    else
+    {
+      NICE::Image img;
+      try
+      {
+        img.read ( filelist[it] );
+      }
+      catch ( ImageException iE )
+      {
+        fprintf ( stderr, "Failed to open image file: %s\n", filelist[it].c_str() );
+        fprintf ( stderr, "%s\n", iE.what() );
+        exit ( -1 );
+      }
+      if ( !isInit )
+      {
+        imgData.reInit ( img.width(),img.height(),filelist.size(),1 );
+        isInit = true;
+      }
+      for ( int y = 0; y < img.height(); y++ )
+      {
+        for ( int x = 0; x < img.width(); x++ )
+        {
+          imgData.set ( x, y, it, img.getPixel ( x,y ), 0 );
+        }
+      }
+    }
+  }
+
+  if ( imagetype == IMAGETYPE_GRAY )
+  {
+    imgData.equalizeHistogram( 0 );
+#ifdef DEBUG_PRINTS
+    for (int z = 0; z < imgData.depth(); z++)
+    {
+      NICE::Image im = imgData.getChannel( z, 0);
+      im.write( filelist[z]+"_eq.pgm");
+    }
+#endif
+  }
+}
+
+void SemanticSegmentation::getProbabilityMap ( const NICE::MultiChannelImage3DT<double> & prob )
+{
+  std::string s;
+
+  for ( int cl = 0; cl < prob.channels(); cl++ )
+    for ( int z = 0; z < prob.depth(); z++ )
+    {
+      NICE::ColorImage img( prob.width(),prob.height() );
+      NICE::ImageT<double> m = prob.getChannelT(z, cl);
+      imageToPseudoColor(m, img);
+
+      std::stringstream out;
+      out << "probmap_s" << z << "_c" << cl << ".ppm";
+      s = out.str();
+      img.write( s );
+
+      //showImage(img, "Probability map");
+      //getchar();
+    }
+}
+
 ///////////////////// INTERFACE PERSISTENT /////////////////////
 // interface specific methods for store and restore
-///////////////////// INTERFACE PERSISTENT ///////////////////// 
+///////////////////// INTERFACE PERSISTENT /////////////////////
 
 void SemanticSegmentation::restore ( std::istream & is, int format )
 {
   //delete everything we knew so far...
   this->clear();
-  
+
   bool b_restoreVerbose ( false );
 #ifdef B_RESTOREVERBOSE
   b_restoreVerbose = true;
-#endif  
-  
+#endif
+
   if ( is.good() )
   {
-    if ( b_restoreVerbose ) 
+    if ( b_restoreVerbose )
       std::cerr << " restore SemanticSegmentation" << std::endl;
-    
+
     std::string tmp;
-    is >> tmp; //class name 
-    
+    is >> tmp; //class name
+
     if ( ! this->isStartTag( tmp, "SemanticSegmentation" ) )
     {
       std::cerr << " WARNING - attempt to restore SemanticSegmentation, but start flag " << tmp << " does not match! Aborting... " << std::endl;
       throw;
-    }   
-    
+    }
+
     is.precision (numeric_limits<double>::digits10 + 1);
-    
+
     bool b_endOfBlock ( false ) ;
-    
+
     while ( !b_endOfBlock )
     {
-      is >> tmp; // start of block 
-      
+      is >> tmp; // start of block
+
       if ( this->isEndTag( tmp, "SemanticSegmentation" ) )
       {
         b_endOfBlock = true;
         continue;
-      }      
-      
+      }
+
       tmp = this->removeStartTag ( tmp );
-      
+
       if ( b_restoreVerbose )
         std::cerr << " currently restore section " << tmp << " in SemanticSegmentation" << std::endl;
-           
+
       if ( tmp.compare("classNames") == 0 )
       {
-	//dirty solution to circumvent the const-flag
+  //dirty solution to circumvent the const-flag
         const_cast<ClassNames*>(this->classNames)->restore ( is, format );
-        
-        is >> tmp; // end of block 
+
+        is >> tmp; // end of block
         tmp = this->removeEndTag ( tmp );
       }
       else if ( tmp.compare("imagetype") == 0 )
       {
         unsigned int ui_imagetyp;
-        is >> ui_imagetyp;        
-        this->imagetype = static_cast<IMAGETYP> ( ui_imagetyp );	
-        
-        is >> tmp; // end of block 
+        is >> ui_imagetyp;
+        this->imagetype = static_cast<IMAGETYP> ( ui_imagetyp );
+
+        is >> tmp; // end of block
         tmp = this->removeEndTag ( tmp );
       }
       else if ( tmp.compare("iterationCountSuffix") == 0 )
       {
         is >> this->iterationCountSuffix;
-        
-        is >> tmp; // end of block 
+
+        is >> tmp; // end of block
         tmp = this->removeEndTag ( tmp );
-      } 
+      }
       else
       {
       std::cerr << "WARNING -- unexpected SemanticSegmentation object -- " << tmp << " -- for restoration... aborting" << std::endl;
@@ -330,38 +476,38 @@ void SemanticSegmentation::restore ( std::istream & is, int format )
     std::cerr << "SemanticSegmentation::restore -- InStream not initialized - restoring not possible!" << std::endl;
     throw;
   }
-  
-  
+
+
  //TODO check whether we also have to do something linke Preprocess::Init ( conf );
 }
 
 void SemanticSegmentation::store ( std::ostream & os, int format ) const
-{ 
+{
   if (os.good())
   {
     // show starting point
-    os << this->createStartTag( "SemanticSegmentation" ) << std::endl;    
-    
+    os << this->createStartTag( "SemanticSegmentation" ) << std::endl;
+
     os.precision (numeric_limits<double>::digits10 + 1);
-    
+
     os << this->createStartTag( "classNames" ) << std::endl;
     this->classNames->store ( os, format );
-    os << this->createEndTag( "classNames" ) << std::endl;    
-    
+    os << this->createEndTag( "classNames" ) << std::endl;
+
     //
-    
+
     os << this->createStartTag( "imagetype" ) << std::endl;
     os << imagetype << std::endl;
-    os << this->createEndTag( "imagetype" ) << std::endl; 
-    
+    os << this->createEndTag( "imagetype" ) << std::endl;
+
     //
-    
+
     os << this->createStartTag( "iterationCountSuffix" ) << std::endl;
     os << iterationCountSuffix << std::endl;
-    os << this->createEndTag( "iterationCountSuffix" ) << std::endl;     
-    
+    os << this->createEndTag( "iterationCountSuffix" ) << std::endl;
+
     // done
-    os << this->createEndTag( "SemanticSegmentation" ) << std::endl;    
+    os << this->createEndTag( "SemanticSegmentation" ) << std::endl;
   }
   else
   {
@@ -373,4 +519,3 @@ void SemanticSegmentation::clear ()
 {
  //TODO
 }
-

semseg/SemanticSegmentation.h

@@ -1,28 +1,31 @@
 /**
 * @file SemanticSegmentation.h
 * @brief abstract interface for semantic segmentation algorithms
-* @author Erik Rodner, Alexander Freytag
-* @date 03/19/2009, latest update: 29-01-2014 (dd-mm-yyyy)
+* @author Erik Rodner, Alexander Freytag, Sven Sickert
+* @date 03/19/2009, latest update: 14-05-2014 (dd-mm-yyyy)
 
 */
 #ifndef SEMANTICSEGMENTATIONINCLUDE
 #define SEMANTICSEGMENTATIONINCLUDE
 
+// standard library includes
+#include <vector>
+
 // nice-core includes
 #include <core/basics/Persistent.h>
+#include "core/image/MultiChannelImage3DT.h"
 
 // nice-vislearning includes
-#include <vislearning/cbaselib/MultiDataset.h>
-#include <vislearning/cbaselib/LocalizationResult.h>
-#include <vislearning/cbaselib/CachedExample.h>
-#include <vislearning/cbaselib/Example.h>
+#include "vislearning/cbaselib/MultiDataset.h"
+#include "vislearning/cbaselib/LocalizationResult.h"
+#include "vislearning/cbaselib/CachedExample.h"
+#include "vislearning/cbaselib/Example.h"
 
 
 #define ROADWORKSADD fthrow(NICE::Exception, "addNewExample(const NICE::Vector & newExample, const int & newClassNo): not yet implemented!");
 #define ROADWORKSADDNOVEL fthrow(NICE::Exception, "addNovelExamples(): not yet implemented!");
 #define ROADWORKSGETNOVEL fthrow(NICE::Exception, "getNovelExamples(): not yet implemented!");
 
-
 namespace OBJREC
 {
 
@@ -31,13 +34,13 @@ class SemanticSegmentation : public NICE::Persistent
 {
 
   protected:
-    
+
     /////////////////////////
     /////////////////////////
     // PROTECTED VARIABLES //
     /////////////////////////
-    ///////////////////////// 
-    
+    /////////////////////////
+
     /** accessible class names and information about
         number of classes etc. */
      const ClassNames * classNames;
@@ -51,27 +54,27 @@ class SemanticSegmentation : public NICE::Persistent
 
     /** whether to load images with color information */
     IMAGETYP imagetype;
-    
+
     int iterationCountSuffix;
-    
+
     /////////////////////////
     /////////////////////////
     //  PROTECTED METHODS  //
     /////////////////////////
-    /////////////////////////    
+    /////////////////////////
 
   public:
-    
+
     ///////////////////// ///////////////////// /////////////////////
     //                   CONSTRUCTORS / DESTRUCTORS
-    ///////////////////// ///////////////////// /////////////////////     
+    ///////////////////// ///////////////////// /////////////////////
 
     /** default constructor
      * @author Alexander Freytag
      * @date 06-02-2014 ( dd-mm-yyy )
     */
-    SemanticSegmentation ( );    
-    
+    SemanticSegmentation ( );
+
     /** simple constructor
         @param conf global settings
         @param classNames this ClassNames object while be stored as a attribute
@@ -81,23 +84,23 @@ class SemanticSegmentation : public NICE::Persistent
 
     /** simple destructor */
     virtual ~SemanticSegmentation();
-    
-    /** 
+
+    /**
     * @brief Setup internal variables and objects used
     * @author Alexander Freytag
     * @param conf Config file to specify variable settings
     * @param s_confSection
-    */    
-    void initFromConfig(const NICE::Config *conf, const std::string & s_confSection = "SemanticSegmentation");  
-    
+    */
+    void initFromConfig(const NICE::Config *conf, const std::string & s_confSection = "SemanticSegmentation");
+
     ///////////////////// ///////////////////// /////////////////////
     //                      SEGMENTATION STUFF
-    ///////////////////// ///////////////////// /////////////////////      
+    ///////////////////// ///////////////////// /////////////////////
 
     /** load img from file call localize(CachedExample *ce) etc. */
     void semanticseg ( const std::string & filename,
                        NICE::Image & segresult,
-                       NICE::MultiChannelImageT<double> & probabilities );    
+                       NICE::MultiChannelImageT<double> & probabilities );
 
     /** this function has to be overloaded by all subclasses
         @param ce image data
@@ -110,10 +113,39 @@ class SemanticSegmentation : public NICE::Persistent
                                NICE::Image & segresult,
                                NICE::MultiChannelImageT<double> & probabilities ) = 0;
 
-                               
+    /**
+     * @brief Classification function (has to be overloaded by all subclasses)
+     * @author Sven Sickert
+     * @param imgData image data
+     * @param segresult result of the semantic segmentation with a label for each pixel
+     * @param probabilities multi-channel image with one channel for each class and
+     *        corresponding probabilities for each pixel
+     * @param filelist filename list of images that represent slices of a stack
+     */
+    virtual void classify ( NICE::MultiChannelImage3DT<double> & imgData,
+                            NICE::MultiChannelImageT<double> & segresult,
+                            NICE::MultiChannelImage3DT<double> & probabilities,
+                            const std::vector<std::string> & filelist ) = 0;
+
+    /** training function (has to be overloaded by all subclasses)
+     * @param md the data set
+     */
+    virtual void train ( const MultiDataset * md ) = 0;
+
+    /**
+     * @brief Load image slices into a MultiChannelImage3DT
+     * @author Sven Sickert
+     * @param filelist filename list of images that represent slices of a stack
+     * @param imgData output
+     */
+    void make3DImage ( const std::vector<std::string> & filelist,
+                       NICE::MultiChannelImage3DT<double> & imgData );
+
+
     ///////////////////// ///////////////////// /////////////////////
     //                      DATA CONVERSION
-    ///////////////////// ///////////////////// /////////////////////                                
+    ///////////////////// ///////////////////// /////////////////////
+
     /**
      * convert different datatypes
      */
@@ -123,41 +155,56 @@ class SemanticSegmentation : public NICE::Persistent
     void convertLSetToExamples ( OBJREC::Examples &examples, OBJREC::LabeledSetVector &lvec, const bool & removeOldDataPointer=false );
     void convertLSetToSparseExamples ( OBJREC::Examples &examples, OBJREC::LabeledSetVector &lvec );
 
-    
 
-    
     ///////////////////// ///////////////////// /////////////////////
     //                       ONLINE LEARNING
-    ///////////////////// ///////////////////// /////////////////////      
-    
+    ///////////////////// ///////////////////// /////////////////////
+
     virtual void addNewExample(const NICE::Vector & newExample, const int & newClassNo)
     {
-      ROADWORKSADD;      
-    };
+      ROADWORKSADD;
+    }
+
     /**
      * @brief Add those examples, which belong to the most novel region seen so far
      *
      * @return void
-     **/    
+     **/
     virtual void addNovelExamples()
     {
-      ROADWORKSADDNOVEL;      
-    };   
-    
+      ROADWORKSADDNOVEL;
+    }
+
     ///////////////////// ///////////////////// /////////////////////
     //                         GET / SET
-    ///////////////////// ///////////////////// /////////////////////      
-    
+    ///////////////////// ///////////////////// /////////////////////
+
     /**
      * @brief Get a pointer to the examples extracted from the most novel region seen so far
      *
      * @return Examples *
-     **/        
+     **/
     virtual const Examples * getNovelExamples() const
     {
-      ROADWORKSGETNOVEL;      
-    };
-    
+      ROADWORKSGETNOVEL;
+    }
+
+    /**
+     * @brief Collect information about the depth of 3d images
+     * @author Sven Sickert
+     * @param Files a labeled set of data
+     * @param depthVec output of depth values
+     * @param run3dseg whether slice counting is necessary or not
+     */
+    void getDepthVector ( const LabeledSet *Files, std::vector<int> & depthVec, const bool run3dseg );
+
+    /**
+     * @brief Save probability maps of all classes to iamge files
+     * @author Sven Sickert
+     * @param prob class probability maps
+     */
+    void getProbabilityMap( const NICE::MultiChannelImage3DT<double> & prob );
+
     /**
      * @author Alexander Freytag
      * @date 06-02-2014 ( dd-mm-yyyy )
@@ -168,27 +215,25 @@ class SemanticSegmentation : public NICE::Persistent
     
     ///////////////////// INTERFACE PERSISTENT /////////////////////
     // interface specific methods for store and restore
-    ///////////////////// INTERFACE PERSISTENT /////////////////////   
-    
-    /** 
+    ///////////////////// INTERFACE PERSISTENT /////////////////////
+
+    /**
      * @brief Load active-segmentation-object from external file (stream)
      * @author Alexander Freytag
-     */     
+     */
     virtual void restore ( std::istream & is, int format = 0 );
-    
-    /** 
+
+    /**
      * @brief Save active-segmentation-object to external file (stream)
      * @author Alexander Freytag
-     */       
+     */
     virtual void store( std::ostream & os, int format = 0 ) const;
-    
-    /** 
+
+    /**
      * @brief Clear active-segmentation-object object
      * @author Alexander Freytag
-     */    
-    virtual void clear ();    
-    
-    
+     */
+    virtual void clear ();
 
 };
 

semseg/libdepend.inc

@@ -1,3 +1,3 @@
 $(call PKG_DEPEND_EXT,OPENMP)
 $(call PKG_DEPEND_INT,segmentation)
-$(call PKG_DEPEND_INT,gp-hik-exp)
+$(call PKG_DEPEND_INT,gp-hik-exp)

semseg/operations/Operations3D.cpp

@@ -0,0 +1,424 @@
+#include "Operations3D.h"
+
+using namespace OBJREC;
+using namespace std;
+using namespace NICE;
+
+Operation3D::Operation3D()
+{
+  init = false;
+  context = false;
+}
+
+void Operation3D::set ( int _x1, int _y1, int _z1, int _x2, int _y2, int _z2, int _channel1, int _channel2, int _featType )
+{
+  x1 = _x1;
+  y1 = _y1;
+	z1 = _z1;
+  x2 = _x2;
+  y2 = _y2;
+	z2 = _z2;
+  channel1 = _channel1;
+  channel2 = _channel2;
+  featType = _featType;
+  init = true;
+}
+
+void Operation3D::setContext ( bool _context )
+{
+  context = _context;
+}
+
+bool Operation3D::getContext()
+{
+  return context;
+}
+
+void Operation3D::setFeatType ( int _featType )
+{
+  featType = _featType;
+}
+
+int Operation3D::getFeatType()
+{
+  return featType;
+}
+
+void Operation3D::getXYZ ( const Features &feats, int &xsize, int &ysize, int &zsize )
+{
+  xsize = feats.feats->width();
+  ysize = feats.feats->height();
+	zsize = feats.feats->depth();
+}
+
+
+void Operation3D::store ( std::ostream & os )
+{
+  os << x1 << " " << x2 << " " << y1 << " " << y2 << " " << z1 << " " << z2 << " " << channel1 << " " << channel2 << " " << featType << std::endl;
+  if ( !init )
+    os << -1 << std::endl;
+  else
+  {
+    if (featType == 3 || featType == 4)
+      os << CONTEXT << std::endl;
+    else
+      os << RAWFEAT << std::endl;
+  }
+}
+
+void Operation3D::restore ( std::istream &is )
+{
+  is >> x1;
+  is >> x2;
+  is >> y1;
+  is >> y2;
+	is >> z1;
+  is >> z2;
+  is >> channel1;
+  is >> channel2;
+  is >> featType;
+
+  int tmp;
+  is >> tmp;
+
+  if ( tmp >= 0 )
+  {
+    if ( tmp == RAWFEAT || tmp == CONTEXT )
+    {
+      init = true;
+    }
+    else
+    {
+      throw ( "no valid ValueAccess" );
+    }
+  }
+  else
+  {
+    init = false;
+  }
+}
+
+std::string Operation3D::writeInfos()
+{
+  std::stringstream ss;
+  ss << " x1: " << x1 << " y1: " << y1 << " z1: " << z1 << " x2: " << x2 << " y2: " << y2 << " z2: " << z2 <<  " c1: " << channel1 << " c2: " << channel2;
+  return ss.str();
+}
+
+
+
+
+//############################## region feature ###############################
+
+
+double RegionFeat::getVal ( const Features &feats, const int &x, const int &y, const int &z )
+{
+  return (*feats.rProbs)[(*feats.feats)(x,y,z,channel1)][channel2];
+}
+
+
+
+
+//############################# two gray values ###############################
+
+
+double Minus::getVal ( const Features &feats, const int &x, const int &y, const int &z )
+{
+  int xsize, ysize, zsize;
+  getXYZ ( feats, xsize, ysize, zsize );
+
+  double v1 = feats.feats->get ( BOUND ( x + x1, 0, xsize - 1 ), BOUND ( y + y1, 0, ysize - 1 ), BOUND ( z + z1, 0, zsize - 1 ), channel1 );
+  double v2 = feats.feats->get ( BOUND ( x + x2, 0, xsize - 1 ), BOUND ( y + y2, 0, ysize - 1 ), BOUND ( z + z2, 0, zsize - 1 ), channel2 );
+
+  return v1-v2;
+}
+
+double MinusAbs::getVal ( const Features &feats, const int &x, const int &y, const int &z )
+{
+  int xsize, ysize, zsize;
+  getXYZ ( feats, xsize, ysize, zsize );
+
+  double v1 = feats.feats->get ( BOUND ( x + x1, 0, xsize - 1 ), BOUND ( y + y1, 0, ysize - 1 ), BOUND ( z + z1, 0, zsize - 1 ), channel1 );
+  double v2 = feats.feats->get ( BOUND ( x + x2, 0, xsize - 1 ), BOUND ( y + y2, 0, ysize - 1 ), BOUND ( z + z2, 0, zsize - 1 ), channel2 );
+
+  return abs(v1-v2);
+}
+
+double Addition::getVal ( const Features &feats, const int &x, const int &y, const int &z )
+{
+  int xsize, ysize, zsize;
+  getXYZ ( feats, xsize, ysize, zsize );
+
+  double v1 = feats.feats->get ( BOUND ( x + x1, 0, xsize - 1 ), BOUND ( y + y1, 0, ysize - 1 ), BOUND ( z + z1, 0, zsize - 1 ), channel1 );
+  double v2 = feats.feats->get ( BOUND ( x + x2, 0, xsize - 1 ), BOUND ( y + y2, 0, ysize - 1 ), BOUND ( z + z2, 0, zsize - 1 ), channel2 );
+
+  return v1+v2;
+}
+
+
+
+
+//############################## one gray value ###############################
+
+
+double Only1::getVal ( const Features &feats, const int &x, const int &y, const int &z )
+{
+  int xsize, ysize, zsize;
+  getXYZ ( feats, xsize, ysize, zsize );
+
+  return feats.feats->get ( BOUND ( x + x1, 0, xsize - 1 ), BOUND ( y + y1, 0, ysize - 1 ), BOUND ( z + z1, 0, zsize - 1 ), channel1 );
+}
+
+
+
+
+//############################ relative positions #############################
+
+
+double RelativeXPosition::getVal ( const Features &feats, const int &x, const int &y, const int &z )
+{
+  int xsize, ysize, zsize;
+  getXYZ ( feats, xsize, ysize, zsize );
+  return ( double ) x / ( double ) xsize;
+}
+
+double RelativeYPosition::getVal ( const Features &feats, const int &x, const int &y, const int &z )
+{
+  int xsize, ysize, zsize;
+  getXYZ ( feats, xsize, ysize, zsize );
+  return ( double ) y / ( double ) ysize;
+}
+
+double RelativeZPosition::getVal ( const Features &feats, const int &x, const int &y, const int &z )
+{
+  int xsize, ysize, zsize;
+  getXYZ ( feats, xsize, ysize, zsize );
+  return ( double ) z / ( double ) zsize;
+}
+
+
+
+
+
+//########################### integral operations #############################
+
+
+
+double GlobalFeats::getVal ( const Features &feats, const int &x, const int &y, const int &z )
+{
+  int xsize, ysize, zsize;
+  getXYZ ( feats, xsize, ysize, zsize );
+
+  return feats.feats->getIntegralValue( 0, 0, 0, xsize - 1, ysize - 1, zsize - 1, channel1 );
+}
+
+void IntegralOps::set ( int _x1, int _y1, int _z1, int _x2, int _y2, int _z2, int _channel1, int _channel2, int _featType )
+{
+  x1 = std::min ( _x1, _x2 );
+  y1 = std::min ( _y1, _y2 );
+  z1 = std::min ( _z1, _z2 );
+  x2 = std::max ( _x1, _x2 );
+  y2 = std::max ( _y1, _y2 );
+  z2 = std::max ( _z1, _z2 );
+  channel1 = _channel1;
+  channel2 = _channel2;
+  featType = _featType;
+  init = true;
+}
+
+double IntegralOps::getVal ( const Features &feats, const int &x, const int &y, const int &z )
+{
+  return feats.feats->getIntegralValue(x + x1, y + y1, z + z1, x + x2, y + y2, z + z2, channel1);
+}
+
+void IntegralCenteredOps::set ( int _x1, int _y1, int _z1, int _x2, int _y2, int _z2, int _channel1, int _channel2, int _featType )
+{
+  x1 = abs ( _x1 );
+  y1 = abs ( _y1 );
+  z1 = abs ( _z1 );
+  x2 = abs ( _x2 );
+  y2 = abs ( _y2 );
+  z2 = abs ( _z2 );
+  channel1 = _channel1;
+  channel2 = _channel2;
+  featType = _featType;
+  init = true;
+}
+
+double IntegralCenteredOps::getVal ( const Features &feats, const int &x, const int &y, const int &z )
+{
+  return feats.feats->getIntegralValue(x - x1, y - y1, z - z1, x + x1, y + y1, z + z1, channel1);
+}
+
+void BiIntegralCenteredOps::set ( int _x1, int _y1, int _z1, int _x2, int _y2, int _z2, int _channel1, int _channel2, int _featType )
+{
+  x1 = std::min ( abs ( _x1 ), abs ( _x2 ) );
+  y1 = std::min ( abs ( _y1 ), abs ( _y2 ) );
+  z1 = std::min ( abs ( _z1 ), abs ( _z2 ) );
+  x2 = std::max ( abs ( _x1 ), abs ( _x2 ) );
+  y2 = std::max ( abs ( _y1 ), abs ( _y2 ) );
+  z2 = std::max ( abs ( _z1 ), abs ( _z2 ) );
+  channel1 = _channel1;
+  channel2 = _channel2;
+  featType = _featType;
+  init = true;
+}
+
+double BiIntegralCenteredOps::getVal ( const Features &feats, const int &x, const int &y, const int &z )
+{
+  return feats.feats->getIntegralValue(x - x1, y - y1, z - z1, x + x1, y + y1, z + z1, channel1 ) - feats.feats->getIntegralValue(x - x2, y - y2, z - z2, x + x2, y + y2, z + z2, channel1);
+}
+
+
+
+
+
+
+//############################ Haar-like features #############################
+
+
+
+double HaarHorizontal::getVal ( const Features &feats, const int &x, const int &y, const int &z )
+{
+  int tlx = x - abs(x1);
+  int tly = y - abs(y1);
+  int tlz = z - abs(z1);
+  int lrx = x + abs(x1);
+  int lry = y + abs(y1);
+  int lrz = z + abs(z1);
+
+  return feats.feats->getIntegralValue(tlx, tly, tlz, lrx, y, lrz, channel1 ) - feats.feats->getIntegralValue(tlx, y, tlz, lrx, lry, lrz, channel1);
+}
+
+double HaarVertical::getVal ( const Features &feats, const int &x, const int &y, const int &z )
+{
+  int tlx = x - abs(x1);
+  int tly = y - abs(y1);
+  int tlz = z - abs(z1);
+  int lrx = x + abs(x1);
+  int lry = y + abs(y1);
+  int lrz = z + abs(z1);
+
+  return feats.feats->getIntegralValue(tlx, tly, tlz, x, lry, lrz, channel1) - feats.feats->getIntegralValue(x, tly, tlz, lrx, lry, lrz, channel1);
+}
+
+double HaarStacked::getVal ( const Features &feats, const int &x, const int &y, const int &z )
+{
+  int tlx = x - abs(x1);
+  int tly = y - abs(y1);
+  int tlz = z - abs(z1);
+  int lrx = x + abs(x1);
+  int lry = y + abs(y1);
+  int lrz = z + abs(z1);
+
+  return feats.feats->getIntegralValue(tlx, tly, tlz, lrx, lry, z, channel1) - feats.feats->getIntegralValue(tlx, tly, z, lrx, lry, lrz, channel1);
+}
+
+double HaarDiagXY::getVal ( const Features &feats, const int &x, const int &y, const int &z )
+{
+  int tlx = x - abs(x1);
+  int tly = y - abs(y1);
+  int tlz = z - abs(z1);
+  int lrx = x + abs(x1);
+  int lry = y + abs(y1);
+  int lrz = z + abs(z1);
+
+  return feats.feats->getIntegralValue(tlx, tly, tlz, x, y, lrz, channel1) + feats.feats->getIntegralValue(x, y, tlz, lrx, lry, lrz, channel1) - feats.feats->getIntegralValue(tlx, y, tlz, x, lry, lrz, channel1) - feats.feats->getIntegralValue(x, tly, tlz, lrx, y, lrz, channel1);
+}
+
+double HaarDiagXZ::getVal ( const Features &feats, const int &x, const int &y, const int &z )
+{
+  int tlx = x - abs(x1);
+  int tly = y - abs(y1);
+  int tlz = z - abs(z1);
+  int lrx = x + abs(x1);
+  int lry = y + abs(y1);
+  int lrz = z + abs(z1);
+
+  return feats.feats->getIntegralValue(tlx, tly, tlz, x, lry, z, channel1) + feats.feats->getIntegralValue(x, tly, z, lrx, lry, lrz, channel1) - feats.feats->getIntegralValue(tlx, tly, z, x, lry, lrz, channel1) - feats.feats->getIntegralValue(x, tly, tlz, lrx, lry, z, channel1);
+}
+
+double HaarDiagYZ::getVal ( const Features &feats, const int &x, const int &y, const int &z )
+{
+  int tlx = x - abs(x1);
+  int tly = y - abs(y1);
+  int tlz = z - abs(z1);
+  int lrx = x + abs(x1);
+  int lry = y + abs(y1);
+  int lrz = z + abs(z1);
+
+  return feats.feats->getIntegralValue(tlx, tly, tlz, lrx, y, z, channel1) + feats.feats->getIntegralValue(tlx, y, z, lrx, lry, lrz, channel1) - feats.feats->getIntegralValue(tlx, tly, z, lrx, y, lrz, channel1) - feats.feats->getIntegralValue(tlx, y, tlz, lrx, lry, z, channel1);
+}
+
+double Haar3Horiz::getVal ( const Features &feats, const int &x, const int &y, const int &z )
+{
+  int tlx = x - abs(x1);
+  int tly = y - abs(y1);
+  int tlz = z - abs(z1);
+  int lrx = x + abs(x1);
+  int lry = y + abs(y1);
+  int lrz = z + abs(z1);
+  int mtly = y - abs(y1);
+  int mlry = y + abs(y1);
+
+  return feats.feats->getIntegralValue(tlx, tly, tlz, lrx, mtly, lrz, channel1) - feats.feats->getIntegralValue(tlx, mtly, tlz, lrx, mlry, lrz, channel1) + feats.feats->getIntegralValue(tlx, mlry, tlz, lrx, lry, lrz, channel1);
+}
+
+double Haar3Vert::getVal ( const Features &feats, const int &x, const int &y, const int &z )
+{
+  int tlx = x - abs(x1);
+  int tly = y - abs(y1);
+  int tlz = z - abs(z1);
+  int lrx = x + abs(x1);
+  int lry = y + abs(y1);
+  int lrz = z + abs(z1);
+  int mtlx = x - abs(x1);
+  int mlrx = x + abs(x1);
+
+  return feats.feats->getIntegralValue(tlx, tly, tlz, mtlx, lry, lrz, channel1) - feats.feats->getIntegralValue(mtlx, tly, tlz, mlrx, lry, lrz, channel1) + feats.feats->getIntegralValue(mlrx, tly, tlz, lrx, lry, lrz, channel1);
+}
+
+double Haar3Stack::getVal ( const Features &feats, const int &x, const int &y, const int &z )
+{
+  int tlx = x - abs(x1);
+  int tly = y - abs(y1);
+  int tlz = z - abs(z1);
+  int lrx = x + abs(x1);
+  int lry = y + abs(y1);
+  int lrz = z + abs(z1);
+  int mtlz = z - abs(z1);
+  int mlrz = z + abs(z1);
+
+  return feats.feats->getIntegralValue(tlx, tly, tlz, lrx, lry, mtlz, channel1) - feats.feats->getIntegralValue(tlx, tly, mtlz, lrx, lry, mlrz, channel1) + feats.feats->getIntegralValue(tlx, tly, mlrz, lrx, lry, lrz, channel1);
+}
+
+
+
+
+//############################### Ray Features ################################
+double RayDiff::getVal( const Features &feats, const int &x, const int &y, const int &z )
+{
+
+  double dist1 = feats.feats->get( x, y, z, feats.feats->channels()-24+z1 );
+  double dist2 = feats.feats->get( x, y, z, feats.feats->channels()-24+z2 );
+
+  if (dist1 != 0)
+    return (dist1-dist2)/dist1;
+  else
+    return 0.0;
+}
+
+double RayDist::getVal( const Features &feats, const int &x, const int &y, const int &z )
+{
+  return feats.feats->get( x, y, z, feats.feats->channels()-24+z1 );
+}
+
+double RayNorm::getVal( const Features &feats, const int &x, const int &y, const int &z )
+{
+  return feats.feats->get( x, y, z, feats.feats->channels()-16+z1 );
+}
+
+double RayOrient::getVal( const Features &feats, const int &x, const int &y, const int &z )
+{
+  return feats.feats->get( x, y, z, feats.feats->channels()-8+z1 );
+}

semseg/operations/Operations3D.h

@@ -0,0 +1,1602 @@
+/**
+* @file Operation3D.h
+* @brief abstract class for any kind of feature extraction from 3d images
+* @author Björn Fröhlich, Sven Sickert
+* @date 24.04.2012
+
+*/
+
+#include "core/image/MultiChannelImageT.h"
+#include "core/image/MultiChannelImage3DT.h"
+
+#define BOUND(x,min,max) (((x)<(min))?(min):((x)>(max)?(max):(x)))
+
+namespace OBJREC {
+
+class Operation3D;
+
+/**
+ * @brief methods for value access
+ **/
+enum ValueTypes
+{
+  RAWFEAT,
+  CONTEXT
+};
+
+/**
+ * @brief feature extraction methods
+ **/
+enum OperationTypes {
+  MINUS,
+  MINUSABS,
+  ADDITION,
+  ONLY1,
+  INTEGRAL,
+  INTEGRALCENT,
+  BIINTEGRALCENT,
+  HAARHORIZ,
+  HAARVERT,
+	HAARSTACK,
+  HAARDIAGXY,
+	HAARDIAGXZ,
+	HAARDIAGYZ,
+  HAAR3HORIZ,
+  HAAR3VERT,
+	HAAR3STACK,
+  RELATIVEXPOSITION,
+  RELATIVEYPOSITION,
+	RELATIVEZPOSITION,
+  GLOBALFEATS,
+  EQUALITY,
+  RAYDIFF,
+  RAYDIST,
+  RAYORIENT,
+  RAYNORM,
+  MINUS_C,
+  MINUSABS_C,
+  ADDITION_C,
+  ONLY1_C,
+  INTEGRAL_C,
+  INTEGRALCENT_C,
+  BIINTEGRALCENT_C,
+  HAARHORIZ_C,
+  HAARVERT_C,
+  HAARSTACK_C,
+  HAARDIAGXY_C,
+  HAARDIAGXZ_C,
+  HAARDIAGYZ_C,
+  HAAR3HORIZ_C,
+  HAAR3VERT_C,
+  HAAR3STACK_C,
+  GLOBALFEATS_C,
+  NBOPERATIONS
+};
+
+/**
+ * @brief node class for context tree
+ **/
+class TreeNode
+{
+
+  public:
+    /** left child node */
+    int left;
+
+    /** right child node */
+    int right;
+
+    /** position of feat for decision */
+    Operation3D *feat;
+
+    /** decision stamp */
+    double decision;
+
+    /** is the node a leaf or not */
+    bool isleaf;
+
+    /** distribution in current node */
+    std::vector<double> dist;
+
+    /** depth of the node in the tree */
+    int depth;
+
+    /** how many pixels are in this node */
+    int featcounter;
+
+    /** unique number */
+    int nodeNumber;
+
+    /** simple constructor */
+    TreeNode() : left ( -1 ), right ( -1 ), feat ( NULL ), decision ( -1.0 ), isleaf ( false ) {}
+
+    /** standard constructor */
+    TreeNode ( int _left, int _right, Operation3D *_feat, double _decision ) : left ( _left ), right ( _right ), feat ( _feat ), decision ( _decision ), isleaf ( false ) {}
+
+    /**
+     * @brief initialize node
+     * @param _depth current depth in tree
+     * @param _numClasses amount of classes (initialize distribution)
+     * @param _nodeNumber unique node number
+     */
+    void init( int _depth, int _numClasses, int _nodeNumber)
+    {
+      depth = _depth;
+      featcounter = 0;
+      dist = std::vector<double> (_numClasses, 0.0);
+      nodeNumber = _nodeNumber;
+    }
+};
+
+
+
+/**
+ * @brief holds all necessary information for feature extraction of 3d images
+ **/
+struct Features {
+  /** simple features like RGB values */
+  NICE::MultiChannelImage3DT<double> *feats;
+
+  /** probabilities for each region */
+  std::vector<std::vector<double> > *rProbs;
+};
+
+
+
+/**
+ * @brief abstract operation class
+ **/
+class Operation3D
+{
+  protected:
+    /** two different points (e.g. for an rectangle or two positions), channels and size  */
+    int x1, y1, z1, x2, y2, z2, channel1, channel2;
+    
+    /** type of feature */
+    int featType;
+    
+    bool init;
+
+    bool context;
+
+  public:
+
+    /** simple constructor */
+    Operation3D();
+
+    /**
+     * @brief set all parameters
+     * @param _x1 position 1
+     * @param _y1 position 1
+		 * @param _z1 position 1
+     * @param _x2 position 2
+     * @param _y2 position 2
+		 * @param _z2 position 2
+     * @param _channel1 channel 1
+     * @param _channel2 channel 2
+     * @param _ftype feature type
+     * @return void nothing
+     **/
+    virtual void set ( int _x1, int _y1, int _z1, int _x2, int _y2, int _z2, int _channel1, int _channel2, int _featType );
+
+    /**
+     * @brief set whether it is a context feature or not
+     * @param _context context boolean
+     * @return void nothing
+     **/
+    void setContext ( bool _context );
+
+    /**
+     * @brief return context information (set by setContext(bool)
+     *
+     * @return bool whether context is used or not
+     **/
+    bool getContext();
+    
+    /**
+     * @brief set type of feature
+     * @param _featType type of feature
+     * @return void nothing
+     **/
+    void setFeatType ( int _featType );
+
+    /**
+     * @brief return context information (set by setContext(bool)
+     *
+     * @return int get feature type
+     **/
+    int getFeatType();   
+
+    /**
+     * @brief abstract interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+		 * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z ) = 0;
+
+    /**
+     * @brief virtual clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone() = 0;
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos();
+
+    /**
+     * @brief exctract current image boarders
+     * @param feats image information
+     * @param xsize width
+     * @param ysize height
+		 * @param zsize depth
+     * @return void
+     **/
+    inline void getXYZ ( const Features &feats, int &xsize, int &ysize, int &zsize );
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps() = 0;
+
+    /**
+     * @brief store all information for current operation in stream
+     * @param os out stream
+     * @return void
+     **/
+    virtual void store ( std::ostream & os );
+
+    /**
+     * @brief restore all information for current operation from stream
+     * @param is in stream
+     * @return void
+     **/
+    virtual void restore ( std::istream & is );
+};
+
+/**
+ * @brief simple equality check ?(A==B)
+ **/
+class RegionFeat: public Operation3D
+{
+  public:
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+		 * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new RegionFeat();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      return "(-)RegionFeat           " + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+      return EQUALITY;
+    }
+};
+
+/**
+ * @brief simple difference operation A-B
+ **/
+class Minus: public Operation3D
+{
+  public:
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+		 * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new Minus();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      std::string out = "Minus";
+
+      if ( context )
+        out = "(C)" + out;
+      else
+        out = "(R)" + out;
+
+      return out + "                " +Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+    if (context)
+      return MINUS_C;
+    else
+      return MINUS;
+    }
+};
+
+/**
+ * @brief simple absolute difference operation |A-B|
+ **/
+class MinusAbs: public Operation3D
+{
+
+  public:
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+		 * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new MinusAbs();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      std::string out = "MinusAbs";
+
+        if ( context )
+          out = "(C)" + out;
+        else
+          out = "(R)" + out;
+
+      return out + "             " + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+      if (context)
+        return MINUSABS_C;
+      else
+        return MINUSABS;
+    }
+};
+
+/**
+ * @brief simple addition operation A+B
+ **/
+class Addition: public Operation3D
+{
+
+  public:
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+		 * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new Addition();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      std::string out = "Addition";
+
+        if ( context )
+          out = "(C)" + out;
+        else
+          out = "(R)" + out;
+
+      return out + "             " + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+      if (context)
+        return ADDITION_C;
+      else
+        return ADDITION;
+    }
+};
+
+/**
+ * @brief simple single element access operation
+ **/
+class Only1: public Operation3D
+{
+  public:
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+		 * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new Only1();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      std::string out = "Only1";
+
+        if ( context )
+          out = "(C)" + out;
+        else
+          out = "(R)" + out;
+
+      return out + "                " + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+      if (context)
+        return ONLY1_C;
+      else
+        return ONLY1;
+    }
+};
+
+/**
+ * @brief get current relative x position
+ **/
+class RelativeXPosition: public Operation3D
+{
+  public:
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+		 * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new RelativeXPosition();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      std::string out = "RelXPos";
+
+        if ( context )
+          out = "(C)" + out;
+        else
+          out = "(R)" + out;
+
+      return out + "              " + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+      return RELATIVEXPOSITION;
+    }
+};
+
+/**
+ * @brief get current relative y position
+ **/
+class RelativeYPosition: public Operation3D
+{
+  public:
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+		 * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new RelativeYPosition();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      std::string out = "RelYPos";
+
+        if ( context )
+          out = "(C)" + out;
+        else
+          out = "(R)" + out;
+
+      return out + "              " + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+      return RELATIVEYPOSITION;
+    }
+};
+
+/**
+ * @brief get current relative z position
+ **/
+class RelativeZPosition: public Operation3D
+{
+  public:
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+		 * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new RelativeZPosition();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      std::string out = "RelZPos";
+
+        if ( context )
+          out = "(C)" + out;
+        else
+          out = "(R)" + out;
+
+      return out + "              " + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+      return RELATIVEZPOSITION;
+    }
+};
+
+/**
+ * @brief uses mean in a window given by (x1,y1,z1) (x2,y2,z2)
+ **/
+class IntegralOps: public Operation3D
+{
+  public:
+    /**
+     * @brief set all parameters
+     * @param _x1 position 1
+     * @param _y1 position 1
+		 * @param _z1 position 1
+     * @param _x2 position 2
+     * @param _y2 position 2
+		 * @param _z2 position 2
+     * @param _channel1 channel 1
+     * @param _channel2 channel 2
+     * @param _ftype feature type
+     * @return void nothing
+     **/
+    virtual void set ( int _x1, int _y1, int _z1, int _x2, int _y2, int _z2, int _channel1, int _channel2, int _ftype );
+
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+		 * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new IntegralOps();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      std::string out = "IntegralOps";
+
+        if ( context )
+          out = "(C)" + out;
+        else
+          out = "(R)" + out;
+
+      return out + "          " + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+      if (context)
+        return INTEGRAL_C;
+      else
+        return INTEGRAL;
+    }
+};
+
+
+/**
+ * @brief like a global bag of words to model the current appearance of classes in an image without local context
+ **/
+class GlobalFeats: public IntegralOps
+{
+
+  public:
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+		 * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new GlobalFeats();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      std::string out = "GlobalFeats";
+
+        if ( context )
+          out = "(C)" + out;
+        else
+          out = "(R)" + out;
+
+      return out + "          " + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+      if ( context )
+        return GLOBALFEATS_C;
+      else
+        return GLOBALFEATS;
+    }
+};
+
+/**
+ * @brief uses mean of Integral image given by x1, y1, z1 with current pixel as center
+ **/
+class IntegralCenteredOps: public IntegralOps
+{
+  public:
+    /**
+     * @brief set all parameters
+     * @param _x1 position 1
+     * @param _y1 position 1
+		 * @param _z1 position 1
+     * @param _x2 position 2
+     * @param _y2 position 2
+		 * @param _z2 position 2
+     * @param _channel1 channel 1
+     * @param _channel2 channel 2
+     * @param _ftype feature type
+     * @return void nothing
+     **/
+    virtual void set ( int _x1, int _y1, int _z1, int _x2, int _y2, int _z2, int _channel1, int _channel2, int _ftype );
+
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+		 * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new IntegralCenteredOps();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      std::string out = "IntegralCenteredOps";
+
+        if ( context )
+          out = "(C)" + out;
+        else
+          out = "(R)" + out;
+
+      return out + "  " + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+      if ( context )
+        return INTEGRALCENT_C;
+      else
+        return INTEGRALCENT;
+    }
+};
+
+/**
+ * @brief uses different of mean of Integral image given by two windows, where (x1,y1,z1) is the width, height & depth of window1 and (x2,y2,z2) of window 2
+ **/
+class BiIntegralCenteredOps: public IntegralCenteredOps
+{
+  public:
+    /**
+     * @brief set all parameters
+     * @param _x1 position 1
+     * @param _y1 position 1
+		 * @param _z1 position 1
+     * @param _x2 position 2
+     * @param _y2 position 2
+		 * @param _z2 position 2
+     * @param _channel1 channel 1
+     * @param _channel2 channel 2
+     * @param _ftype feature type
+     * @return void nothing
+     **/
+    virtual void set ( int _x1, int _y1, int _z1, int _x2, int _y2, int _z2, int _channel1, int _channel2, int ftype );
+
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+		 * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new BiIntegralCenteredOps();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      std::string out = "BiIntegralCenteredOps";
+
+        if ( context )
+          out = "(C)" + out;
+        else
+          out = "(R)" + out;
+
+      return out + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+      if ( context )
+        return BIINTEGRALCENT_C;
+      else
+        return BIINTEGRALCENT;
+    }
+};
+
+/**
+ * @brief horizontal Haar features
+ * ++
+ * --
+ **/
+class HaarHorizontal: public IntegralCenteredOps
+{
+  public:
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+		 * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new HaarHorizontal();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      std::string out = "HaarHorizontal";
+
+        if ( context )
+          out = "(C)" + out;
+        else
+          out = "(R)" + out;
+
+      return out + "       " + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+      if ( context )
+        return HAARHORIZ_C;
+      else
+        return HAARHORIZ;
+    }
+};
+
+/**
+ * @brief vertical Haar features
+ * +-
+ * +-
+ **/
+class HaarVertical: public IntegralCenteredOps
+{
+  public:
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+		 * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new HaarVertical();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      std::string out = "HaarVertical";
+
+        if ( context )
+          out = "(C)" + out;
+        else
+          out = "(R)" + out;
+
+      return out + "         " + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+    if ( context )
+      return HAARVERT_C;
+    else
+      return HAARVERT;
+    }
+};
+
+/**
+ * @brief stacked (depth) Haar features
+ * +-
+ * +-
+ **/
+class HaarStacked: public IntegralCenteredOps
+{
+  public:
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+		 * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new HaarStacked();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      std::string out = "HaarStacked";
+
+        if ( context )
+          out = "(C)" + out;
+        else
+          out = "(R)" + out;
+
+      return out + "          " + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+    if ( context )
+      return HAARSTACK_C;
+    else
+      return HAARSTACK;
+    }
+};
+
+/**
+ * @brief x-y diagonal Haar features
+ * +-
+ * -+
+ **/
+class HaarDiagXY: public IntegralCenteredOps
+{
+  public:
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+		 * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new HaarDiagXY();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      std::string out = "HaarDiagXY";
+
+        if ( context )
+          out = "(C)" + out;
+        else
+          out = "(R)" + out;
+
+      return out + "           " + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+    if ( context )
+      return HAARDIAGXY_C;
+    else
+      return HAARDIAGXY;
+    }
+};
+
+/**
+ * @brief x-z diagonal Haar features
+ * +-
+ * -+
+ **/
+class HaarDiagXZ: public IntegralCenteredOps
+{
+  public:
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+		 * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new HaarDiagXZ();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      std::string out = "HaarDiagXZ";
+
+        if ( context )
+          out = "(C)" + out;
+        else
+          out = "(R)" + out;
+
+      return out + "           " + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+    if ( context )
+      return HAARDIAGXZ_C;
+    else
+      return HAARDIAGXZ;
+    }
+};
+
+/**
+ * @brief y-z diagonal Haar features
+ * +-
+ * -+
+ **/
+class HaarDiagYZ: public IntegralCenteredOps
+{
+  public:
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+		 * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new HaarDiagYZ();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      std::string out = "HaarDiagYZ";
+
+        if ( context )
+          out = "(C)" + out;
+        else
+          out = "(R)" + out;
+
+      return out + "           " + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+    if ( context )
+      return HAARDIAGYZ_C;
+    else
+      return HAARDIAGYZ;
+    }
+};
+
+/**
+ * @brief horizontal Haar features
+ * +++
+ * ---
+ * +++
+ */
+
+class Haar3Horiz: public BiIntegralCenteredOps
+{
+  public:
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+		 * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new Haar3Horiz();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      std::string out = "Haar3Horiz";
+
+        if ( context )
+          out = "(C)" + out;
+        else
+          out = "(R)" + out;
+
+      return out + "           " + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+    if ( context )
+      return HAAR3HORIZ_C;
+    else
+      return HAAR3HORIZ;
+    }
+};
+
+/**
+ * @brief vertical Haar features
+ * +-+
+ * +-+
+ * +-+
+ */
+class Haar3Vert: public BiIntegralCenteredOps
+{
+  public:
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+		 * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new Haar3Vert();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      std::string out = "Haar3Vert";
+
+        if ( context )
+          out = "(C)" + out;
+        else
+          out = "(R)" + out;
+
+      return out + "            " + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+    if ( context )
+      return HAAR3VERT_C;
+    else
+      return HAAR3VERT;
+    }
+};
+
+/**
+ * @brief stacked Haar features
+ * +-+
+ * +-+
+ * +-+
+ */
+class Haar3Stack: public BiIntegralCenteredOps
+{
+  public:
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+		 * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new Haar3Stack();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      std::string out = "Haar3Stack";
+
+        if ( context )
+          out = "(C)" + out;
+        else
+          out = "(R)" + out;
+
+      return out + "           " + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+    if ( context )
+      return HAAR3STACK_C;
+    else
+      return HAAR3STACK;
+    }
+};
+
+/**
+ * @brief Ray features
+ */
+/**
+ * @brief Ray Distance-Difference features
+ */
+class RayDiff: public Operation3D
+{
+  public:
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+     * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new RayDiff();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      return "(-)RayDiff              " + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+      return RAYDIFF;
+    }
+};
+
+/**
+ * @brief Ray Distance features
+ */
+class RayDist: public Operation3D
+{
+  public:
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+     * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new RayDist();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      return "(-)RayDist              " + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+      return RAYDIST;
+    }
+};
+
+/**
+ * @brief Ray Orientation features
+ */
+class RayOrient: public Operation3D
+{
+  public:
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+     * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new RayOrient();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      return "(-)RayOrient            " + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+      return RAYORIENT;
+    }
+};
+
+/**
+ * @brief Ray Norm features
+ */
+class RayNorm: public Operation3D
+{
+  public:
+    /**
+     * @brief interface for feature computation
+     * @param feats features
+     * @param x current x position
+     * @param y current y position
+     * @param z current z position
+     * @return double distance
+     **/
+    virtual double getVal ( const Features &feats, const int &x, const int &y, const int &z );
+
+    /**
+     * @brief clone operation instead of copy constructor (copy constructor does not work)
+     **/
+    virtual Operation3D* clone()
+    {
+      return new RayNorm();
+    }
+
+    /**
+     * @brief print some infos about operation extraction type
+     * @return string feature type
+     **/
+    virtual std::string writeInfos()
+    {
+      return "(-)RayNorm              " + Operation3D::writeInfos();
+    }
+
+    /**
+     * @brief return operation type (for store and restore)
+     * @return OperationTypes
+     **/
+    virtual OperationTypes getOps()
+    {
+      return RAYNORM;
+    }
+};
+
+} //end namespace