NICE_SemSeg/semseg/SemSegContextTree.cpp

#include "SemSegContextTree.h"
#include "vislearning/baselib/Globals.h"
#include "vislearning/baselib/ProgressBar.h"
#include "core/basics/StringTools.h"

#include "core/imagedisplay/ImageDisplay.h"

#include "vislearning/cbaselib/CachedExample.h"
#include "vislearning/cbaselib/PascalResults.h"
#include "vislearning/baselib/ColorSpace.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/image/FilterT.h"

#include <omp.h>
#include <iostream>

//#define DEBUG

using namespace OBJREC;
using namespace std;
using namespace NICE;

SemSegContextTree::SemSegContextTree ( const Config *conf, const MultiDataset *md )
    : SemanticSegmentation ( conf, & ( md->getClassNames ( "train" ) ) )
{
  this->conf = conf;
  string section = "SSContextTree";
  lfcw = new LFColorWeijer ( conf );
  firstiteration = true;

  maxSamples = conf->gI ( section, "max_samples", 2000 );

  minFeats = conf->gI ( section, "min_feats", 50 );

  maxDepth = conf->gI ( section, "max_depth", 10 );

  windowSize = conf->gI ( section, "window_size", 16 );

  featsPerSplit = conf->gI ( section, "feats_per_split", 200 );

  useShannonEntropy = conf->gB ( section, "use_shannon_entropy", true );

  nbTrees = conf->gI ( section, "amount_trees", 1 );

  string segmentationtype = conf->gS ( section, "segmentation_type", "slic" );

  useCategorization = conf->gB ( section, "use_categorization", false );

  cndir = conf->gS ( "SSContextTree", "cndir", "" );

  if ( useCategorization && cndir == "" )
  {
    fasthik = new GPHIKClassifierNICE ( conf );
  }
  else
  {
    fasthik = NULL;
  }

  randomTests = conf->gI ( section, "random_tests", 10 );

  bool saveLoadData = conf->gB ( "debug", "save_load_data", false );
  string fileLocation = conf->gS ( "debug", "datafile", "tmp.txt" );

  pixelWiseLabeling = false;

  useRegionFeature = conf->gB ( section, "use_region_feat", false );

  if ( segmentationtype == "meanshift" )
    segmentation = new RSMeanShift ( conf );
  else if ( segmentationtype == "none" )
  {
    segmentation = NULL;
    pixelWiseLabeling = true;
    useRegionFeature = false;
  }
  else if ( segmentationtype == "felzenszwalb" )
    segmentation = new RSGraphBased ( conf );
  else if ( segmentationtype == "slic" )
    segmentation = new RSSlic ( conf );
  else
    throw ( "no valid segmenation_type\n please choose between none, meanshift, slic and felzenszwalb\n" );

  ftypes = conf->gI ( section, "features", 100 );;

  string featsec = "Features";

  vector<Operation*> tops;

  if ( conf->gB ( featsec, "minus", true ) )
    tops.push_back ( new Minus() );
  if ( conf->gB ( featsec, "minus_abs", true ) )
    tops.push_back ( new MinusAbs() );
  if ( conf->gB ( featsec, "addition", true ) )
    tops.push_back ( new Addition() );
  if ( conf->gB ( featsec, "only1", true ) )
    tops.push_back ( new Only1() );
  if ( conf->gB ( featsec, "rel_x", true ) )
    tops.push_back ( new RelativeXPosition() );
  if ( conf->gB ( featsec, "rel_y", true ) )
    tops.push_back ( new RelativeYPosition() );
  if ( conf->gB ( featsec, "rel_z", true ) )
    tops.push_back ( new RelativeZPosition() );

  ops.push_back ( tops );

  tops.clear();
  tops.push_back ( new RegionFeat() );
  ops.push_back ( tops );

  tops.clear();
  if ( conf->gB ( featsec, "int", true ) )
    tops.push_back ( new IntegralOps() );
  if ( conf->gB ( featsec, "bi_int_cent", true ) )
    tops.push_back ( new BiIntegralCenteredOps() );
  if ( conf->gB ( featsec, "int_cent", true ) )
    tops.push_back ( new IntegralCenteredOps() );
  if ( conf->gB ( featsec, "haar_horz", true ) )
    tops.push_back ( new HaarHorizontal() );
  if ( conf->gB ( featsec, "haar_vert", true ) )
    tops.push_back ( new HaarVertical() );
  if ( conf->gB ( featsec, "haar_stack", true ) )
    tops.push_back ( new HaarStacked() );
  if ( conf->gB ( featsec, "haar_diagxy", true ) )
    tops.push_back ( new HaarDiagXY() );
  if ( conf->gB ( featsec, "haar_diagxz", true ) )
    tops.push_back ( new HaarDiagXZ() );
  if ( conf->gB ( featsec, "haar_diagyz", true ) )
    tops.push_back ( new HaarDiagYZ() );
  if ( conf->gB ( featsec, "haar3_horz", true ) )
    tops.push_back ( new Haar3Horiz() );
  if ( conf->gB ( featsec, "haar3_vert", true ) )
    tops.push_back ( new Haar3Vert() );
  if ( conf->gB ( featsec, "haar3_stack", true ) )
    tops.push_back ( new Haar3Stack() );

  ops.push_back ( tops );
  ops.push_back ( tops );

  tops.clear();
  if ( conf->gB ( featsec, "minus", true ) )
    tops.push_back ( new Minus() );
  if ( conf->gB ( featsec, "minus_abs", true ) )
    tops.push_back ( new MinusAbs() );
  if ( conf->gB ( featsec, "addition", true ) )
    tops.push_back ( new Addition() );
  if ( conf->gB ( featsec, "only1", true ) )
    tops.push_back ( new Only1() );
  if ( conf->gB ( featsec, "rel_x", true ) )
    tops.push_back ( new RelativeXPosition() );
  if ( conf->gB ( featsec, "rel_y", true ) )
    tops.push_back ( new RelativeYPosition() );
  if ( conf->gB ( featsec, "rel_z", true ) )
    tops.push_back ( new RelativeZPosition() );

  ops.push_back ( tops );

  useGradient = conf->gB ( featsec, "use_gradient", true );

  useWeijer = conf->gB ( featsec, "use_weijer", true );

  useGaussian = conf->gB ( featsec, "use_diff_gaussian", false );
  
  useVariance = conf->gB ( featsec, "use_variance", false );

  // geometric features of hoiem
  useHoiemFeatures = conf->gB ( featsec, "use_hoiem_features", false );
  if ( useHoiemFeatures )
  {
    hoiemDirectory = conf->gS ( featsec, "hoiem_directory" );
  }

  opOverview = vector<int> ( NBOPERATIONS, 0 );
  contextOverview = vector<vector<double> > ( maxDepth, vector<double> ( 2, 0.0 ) );

  calcVal.push_back ( new MCImageAccess() );
  calcVal.push_back ( new MCImageAccess() );
  calcVal.push_back ( new MCImageAccess() );
  calcVal.push_back ( new MCImageAccess() );
  calcVal.push_back ( new ClassificationResultAccess() );


  classnames = md->getClassNames ( "train" );

  ///////////////////////////////////
  // Train Segmentation Context Trees
  ///////////////////////////////////

  if ( saveLoadData )
  {
    if ( FileMgt::fileExists ( fileLocation ) )
      read ( fileLocation );
    else
    {
      train ( md );
      write ( fileLocation );
    }
  }
  else
  {
    train ( md );
  }
}

SemSegContextTree::~SemSegContextTree()
{
}

double SemSegContextTree::getBestSplit ( std::vector<NICE::MultiChannelImage3DT<double> > &feats, std::vector<NICE::MultiChannelImage3DT<unsigned short int> > &currentfeats, const std::vector<NICE::MultiChannelImageT<int> > &labels, int node, Operation *&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;

  set<vector<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 ) currentfeats[iCounter].width();
    int ysize = ( int ) currentfeats[iCounter].height();
    int zsize = ( int ) currentfeats[iCounter].depth();

    for ( int x = 0; x < xsize; x++ )
    {
      for ( int y = 0; y < ysize; y++ )
      {
        for ( int z = 0; z < zsize; z++ )
        {
          if ( currentfeats[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]] )
            {
              vector<int> tmp ( 4, 0 );
              tmp[0] = iCounter;
              tmp[1] = x;
              tmp[2] = y;
              tmp[3] = z;
              featcounter++;
              selFeats.insert ( tmp );
              e[cn]++;
            }
          }
        }
      }
    }
  }

  map<int, int>::iterator mapit;

  // global entropy
  double globent = 0.0;

  for ( 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;
  }

  // vector of all possible features
  std::vector<Operation*> featsel;

  for ( int i = 0; i < featsPerSplit; i++ )
  {
    int x1, x2, y1, y2, z1, z2;
    int ft = ( int ) ( ( double ) rand() / ( double ) RAND_MAX * ( double ) ftypes );

    int tmpws = windowSize;

    if ( firstiteration )
      ft = 0;

    if ( channelsPerType[ft].size() == 0 )
    {
      ft = 0;
    }

    if ( ft > 1 )
    {
      //use larger window size for context features
      tmpws *= 2;
    }

    if ( ft == 1 )
    {
      if ( depth < 8 )
      {
        ft = 0;
      }
    }

    // random value range between (-window-size/2) and (window-size/2)
    double z_ratio = conf->gB ( "SSContextTree", "z_ratio", 1.0 );
    int tmp_z =  ( int ) floor( (tmpws * z_ratio) + 0.5 );
    x1 = ( int ) ( ( double ) rand() / ( double ) RAND_MAX * ( double ) tmpws ) - tmpws / 2;
    x2 = ( int ) ( ( double ) rand() / ( double ) RAND_MAX * ( double ) tmpws ) - tmpws / 2;
    y1 = ( int ) ( ( double ) rand() / ( double ) RAND_MAX * ( double ) tmpws ) - tmpws / 2;
    y2 = ( int ) ( ( double ) rand() / ( double ) RAND_MAX * ( double ) tmpws ) - tmpws / 2;
    z1 = ( int ) ( ( double ) rand() / ( double ) RAND_MAX * ( double ) tmp_z ) - tmp_z / 2;
    z2 = ( int ) ( ( double ) rand() / ( double ) RAND_MAX * ( double ) tmp_z ) - tmp_z / 2;

    int f1 = ( int ) ( ( double ) rand() / ( double ) RAND_MAX * ( double ) channelsPerType[ft].size() );
    int f2 = f1;
    if ( ( double ) rand() / ( double ) RAND_MAX > 0.5 )
      f2 = ( int ) ( ( double ) rand() / ( double ) RAND_MAX * ( double ) channelsPerType[ft].size() );
    int o = ( int ) ( ( double ) rand() / ( double ) RAND_MAX * ( double ) ops[ft].size() );

    f1 = channelsPerType[ft][f1];
    f2 = channelsPerType[ft][f2];
    if ( ft == 1 )
    {
      int classes = ( int ) regionProbs[0][0].size();
      f2 = ( int ) ( ( double ) rand() / ( double ) RAND_MAX * ( double ) classes );
    }

    // only depth-values in front of the current pixel are allowed
    bool z_negative_only = conf->gB ( "SSContextTree", "z_negative_only", false );
    if (z_negative_only)
    {
      z1 = -abs(z1);
      z2 = -abs(z2);
    }

    Operation *op = ops[ft][o]->clone();
    op->set ( x1, y1, z1, x2, y2, z2, f1, f2, calcVal[ft] );
    op->setFeatType ( ft );

    if ( ft == 3 || ft == 4 )
      op->setContext ( true );
    else
      op->setContext ( false );

    featsel.push_back ( op );
  }

  for ( int f = 0; f < featsPerSplit; f++ )
  {
    double l_bestig = -numeric_limits< double >::max();
    double l_splitval = -1.0;
    set<vector<int> >::iterator it;
    vector<double> vals;

    double maxval = -numeric_limits<double>::max();
    double minval = numeric_limits<double>::max();
    for ( it = selFeats.begin() ; it != selFeats.end(); it++ )
    {
      Features feat;
      feat.feats = &feats[ ( *it ) [0]];
      feat.cfeats = &currentfeats[ ( *it ) [0]];
      feat.cTree = tree;
      feat.tree = &forest[tree];

      assert ( forest.size() > ( uint ) tree );
      assert ( forest[tree][0].dist.size() > 0 );

      feat.rProbs = &regionProbs[ ( *it ) [0]];

      double val = featsel[f]->getVal ( feat, ( *it ) [1], ( *it ) [2], ( *it ) [3] );
      if ( !isfinite ( val ) )
      {
        //cerr << "feat " << feat.feats->width() << " " << feat.feats->height() << " " << feat.feats->depth() << endl;
        //cerr << "non finite value " << val << " for " << featsel[f]->writeInfos() <<  endl << (*it) [1] << " " <<  (*it) [2] << " " << (*it) [3] << endl;
        val = 0.0;
      }
      vals.push_back ( val );
      maxval = std::max ( val, maxval );
      minval = std::min ( val, minval );
    }

    if ( minval == maxval )
      continue;

    double scale = maxval - minval;
    vector<double> splits;

    for ( int r = 0; r < randomTests; r++ )
    {
      splits.push_back ( ( ( double ) rand() / ( double ) RAND_MAX*scale ) + minval );
    }

    for ( int run = 0 ; run < randomTests; run++ )
    {
      set<vector<int> >::iterator it2;
      double val = splits[run];

      map<int, int> eL, eR;
      int counterL = 0, counterR = 0;
      int counter2 = 0;

      for ( it2 = selFeats.begin() ; it2 != selFeats.end(); it2++, counter2++ )
      {
        int cn = labels[ ( *it2 ) [0]].get ( ( *it2 ) [1], ( *it2 ) [2], ( *it2 ) [3] );
        //cout << "vals[counter2] " << vals[counter2] << " val: " <<  val << endl;

        if ( vals[counter2] < val )
        {
          //left entropie:
          eL[cn] = eL[cn] + 1;
          counterL++;
        }
        else
        {
          //right entropie:
          eR[cn] = eR[cn] + 1;
          counterR++;
        }
      }

      double leftent = 0.0;

      for ( mapit = eL.begin() ; mapit != eL.end(); mapit++ )
      {
        double p = ( double ) ( *mapit ).second / ( double ) counterL;
        leftent -= p * log2 ( p );
      }

      double rightent = 0.0;

      for ( 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 );

      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 = val;
      }
    }

    if ( l_bestig > bestig )
    {
      bestig = l_bestig;
      splitop = featsel[f];
      splitval = l_splitval;
    }
  }
  //FIXME: delete all features!
  /*for(int i = 0; i < featsPerSplit; i++)
  {
   if(featsel[i] != splitop)
    delete featsel[i];
  }*/


#ifdef DEBUG
  //cout << "globent: " << globent <<  " bestig " << bestig << " splitval: " << splitval << endl;
#endif
  return bestig;
}

inline double SemSegContextTree::getMeanProb ( const int &x, const int &y, const int &z, const int &channel, const MultiChannelImage3DT<unsigned short int> &currentfeats )
{
  double val = 0.0;

  for ( int tree = 0; tree < nbTrees; tree++ )
  {
    val += forest[tree][currentfeats.get ( x,y,z,tree ) ].dist[channel];
  }

  return val / ( double ) nbTrees;
}

void SemSegContextTree::computeIntegralImage ( const NICE::MultiChannelImage3DT<unsigned short int> &currentfeats, NICE::MultiChannelImage3DT<double> &feats, int firstChannel )
{
  int xsize = feats.width();
  int ysize = feats.height();
  int zsize = feats.depth();

  if ( firstiteration )
  {
#pragma omp parallel for
    for ( int it = 0; it < ( int ) integralMap.size(); it++ )
    {
      int corg = integralMap[it].first;
      int cint = integralMap[it].second;

      for ( int z = 0; z < zsize; z++ )
      {
        for ( int y = 0; y < ysize; y++ )
        {
          for ( int x = 0; x < xsize; x++ )
          {
            feats ( x, y, z, cint ) = feats ( x, y, z, corg );
          }
        }
      }
      feats.calcIntegral ( cint );
    }
  }

  int channels = ( int ) forest[0][0].dist.size();

#pragma omp parallel for
  for ( int c = 0; c < channels; c++ )
  {

    feats ( 0, 0, 0, firstChannel + c ) = getMeanProb ( 0, 0, 0, c, currentfeats );

    //first column
    for ( int y = 1; y < ysize; y++ )
    {
      feats ( 0, y, 0, firstChannel + c ) = getMeanProb ( 0, y, 0, c, currentfeats )
                                            + feats ( 0, y - 1, 0, firstChannel + c );
    }

    //first row
    for ( int x = 1; x < xsize; x++ )
    {
      feats ( x, 0, 0, firstChannel + c ) = getMeanProb ( x, 0, 0, c, currentfeats )
                                            + feats ( x - 1, 0, 0, firstChannel + c );
    }

    //first stack
    for ( int z = 1; z < zsize; z++ )
    {
      feats ( 0, 0, z, firstChannel + c ) = getMeanProb ( 0, 0, z, c, currentfeats )
                                            + feats ( 0, 0, z - 1, firstChannel + c );
    }

    //x-y plane
    for ( int y = 1; y < ysize; y++ )
    {
      for ( int x = 1; x < xsize; x++ )
      {
        feats ( x, y, 0, firstChannel + c ) = getMeanProb ( x, y, 0, c, currentfeats )
                                              + feats ( x, y - 1, 0, firstChannel + c )
                                              + feats ( x - 1, y, 0, firstChannel + c )
                                              - feats ( x - 1, y - 1, 0, firstChannel + c );
      }
    }

    //y-z plane
    for ( int z = 1; z < zsize; z++ )
    {
      for ( int y = 1; y < ysize; y++ )
      {
        feats ( 0, y, z, firstChannel + c ) = getMeanProb ( 0, y, z, c, currentfeats )
                                              + feats ( 0, y - 1, z, firstChannel + c )
                                              + feats ( 0, y, z - 1, firstChannel + c )
                                              - feats ( 0, y - 1, z - 1, firstChannel + c );
      }
    }

    //x-z plane
    for ( int z = 1; z < zsize; z++ )
    {
      for ( int x = 1; x < xsize; x++ )
      {
        feats ( x, 0, z, firstChannel + c ) = getMeanProb ( x, 0, z, c, currentfeats )
                                              + feats ( x - 1, 0, z, firstChannel + c )
                                              + feats ( x, 0, z - 1, firstChannel + c )
                                              - feats ( x - 1, 0, z - 1, firstChannel + c );
      }
    }

    //rest
    for ( int z = 1; z < zsize; z++ )
    {
      for ( int y = 1; y < ysize; y++ )
      {
        for ( int x = 1; x < xsize; x++ )
        {
          feats ( x, y, z, firstChannel + c ) = getMeanProb ( x, y, z, c, currentfeats )
                                                + feats ( x - 1, y, z, firstChannel + c )
                                                + feats ( x, y - 1, z, firstChannel + c )
                                                + feats ( x, y, z - 1, firstChannel + c )
                                                + feats ( x - 1, y - 1, z - 1, firstChannel + c )
                                                - feats ( x - 1, y - 1, z, firstChannel + c )
                                                - feats ( x - 1, y, z - 1, firstChannel + c )
                                                - feats ( x, y - 1, z - 1, firstChannel + c );
        }
      }
    }
  }
}

inline double computeWeight ( const double &d, const double &dim )
{
  return 1.0 / ( pow ( 2, ( double ) ( dim - d + 1 ) ) );
}



void SemSegContextTree::train ( const MultiDataset *md )
{
  int shortsize = numeric_limits<short>::max();

  Timer timer;
  timer.start();
  const LabeledSet train = * ( *md ) ["train"];
  const LabeledSet *trainp = &train;

  vector<int> zsizeVec;
  getDepthVector ( &train, zsizeVec );
  bool run_3dseg = conf->gB ( "debug", "run_3dseg", true );

  ProgressBar pb ( "compute feats" );
  pb.show();

  //TODO: Speichefresser!, lohnt sich sparse?
  vector<MultiChannelImage3DT<double> > allfeats;
  vector<MultiChannelImage3DT<unsigned short int> > currentfeats;
  vector<MultiChannelImageT<int> > labels;

  vector<SparseVector*> globalCategorFeats;
  vector<map<int,int> > classesPerImage;

  std::string forbidden_classes_s = conf->gS ( "analysis", "donttrain", "" );

  vector<vector<vector<double> > > regionProbs;
  vector<vector<int> > rSize;
  vector<int> amountRegionpI;

  if ( forbidden_classes_s == "" )
  {
    forbidden_classes_s = conf->gS ( "analysis", "forbidden_classes", "" );
  }

  classnames.getSelection ( forbidden_classes_s, forbidden_classes );

  int imgcounter = 0;

  int amountPixels = 0;

  ////////////////////////////////////////////////////
  //define which featurextraction methods should be used for each channel
  if ( imagetype == IMAGETYPE_RGB )
  {
    rawChannels = 3;
  }
  else
  {
    rawChannels = 1;
  }

  // how many channels without integral image
  int shift = 0;

  if ( useGradient )
    rawChannels *= 2;

  if ( useWeijer )
    rawChannels += 11;

  if ( useHoiemFeatures )
    rawChannels += 8;

  if ( useGaussian )
    rawChannels += 1;
  
  if ( useVariance )
    rawChannels += 1; 

  // gray value images
  for ( int i = 0; i < rawChannels; i++ )
  {
    channelType.push_back ( 0 );
  }

  // regions
  if ( useRegionFeature )
  {
    channelType.push_back ( 1 );
    shift++;
  }

///////////////////////////// read input data /////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

  int depthCount = 0;
  vector< string > filelist;
  NICE::MultiChannelImageT<uchar> pixelLabels;

  LOOP_ALL_S ( *trainp )
  {
    EACH_INFO ( classno, info );
    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 ( run_3dseg )
    {
      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 );

    currentfeats.push_back ( MultiChannelImage3DT<unsigned short int> ( xsize, ysize, zsize, nbTrees ) );
    currentfeats[imgcounter].setAll ( 0 );

    //TODO: resize image?!
    MultiChannelImage3DT<double> feats;
    allfeats.push_back ( feats );

    int amountRegions;
    // read image and do some simple transformations
    extractBasicFeatures ( allfeats[imgcounter], imgData, filelist, amountRegions );

    if ( useRegionFeature )
    {
      amountRegionpI.push_back ( amountRegions );
      rSize.push_back ( 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[imgcounter][allfeats[imgcounter] ( x, y, z, rawChannels ) ]++;
          }
        }
      }
    }

    for ( int x = 0; x < xsize; x++ )
    {
      for ( int y = 0; y < ysize; y++ )
      {
        for ( int z = 0; z < zsize; z++ )
        {
          if ( run_3dseg )
            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 )
    {
      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 ( run_3dseg )
              classno = pixelLabels ( x, y, ( uint ) z );
            else
              classno = pL.getPixelQuick ( x,y );

            if ( forbidden_classes.find ( classno ) != forbidden_classes.end() )
              continue;

            classesPerImage[imgcounter][classno] = 1;
          }
        }
      }
    }

    pb.update ( trainp->count() );
    filelist.clear();
    pixelLabels.reInit ( 0,0,0 );
    depthCount = 0;
    imgcounter++;
  }

  pb.hide();

  map<int, int>::iterator mapit;
  int classes = 0;

  for ( mapit = labelcounter.begin(); mapit != labelcounter.end(); mapit++ )
  {
    labelmap[mapit->first] = classes;
    labelmapback[classes] = mapit->first;
    classes++;
  }

//////////////////////////// channel configuration ////////////////////////////
///////////////////////////////////////////////////////////////////////////////

  for ( int i = 0; i < rawChannels; i++ )
  {
    channelType.push_back ( 2 );
  }

  // integral images
  for ( int i = 0; i < classes; i++ )
  {
    channelType.push_back ( 3 );
  }

  integralMap.clear();
  for ( int i = 0; i < rawChannels; i++ )
  {
    integralMap.push_back ( pair<int, int> ( i, i + rawChannels + shift ) );
  }

  int amountTypes = 5;

  channelsPerType = vector<vector<int> > ( amountTypes, vector<int>() );

  for ( int i = 0; i < ( int ) channelType.size(); i++ )
  {
    channelsPerType[channelType[i]].push_back ( i );
  }

  for ( int i = 0; i < classes; i++ )
  {
    channelsPerType[channelsPerType.size()-1].push_back ( i );
  }

  ftypes = std::min ( amountTypes, ftypes );

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

  if ( useRegionFeature )
  {
    for ( int a = 0; a < ( int ) amountRegionpI.size(); a++ )
    {
      regionProbs.push_back ( vector<vector<double> > ( amountRegionpI[a], vector<double> ( classes, 0.0 ) ) );
    }
  }

  //balancing
  int featcounter = 0;

  a = vector<double> ( classes, 0.0 );

  for ( int iCounter = 0; iCounter < imgcounter; iCounter++ )
  {
    int xsize = ( int ) currentfeats[iCounter].width();
    int ysize = ( int ) currentfeats[iCounter].height();
    int zsize = ( int ) currentfeats[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 DEBUG
  for ( int i = 0; i < ( int ) a.size(); i++ )
  {
    cout << "a[" << i << "]: " << a[i] << endl;
  }

  cout << "a.size: " << a.size() << endl;

#endif

  depth = 0;

  uniquenumber = 0;

  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 allleaf = false;
  //int baseFeatSize = allfeats[0].size();

  timer.stop();
  cerr << "preprocessing finished in: " << timer.getLastAbsolute() << " seconds" << endl;
  timer.start();

  while ( !allleaf && depth < maxDepth )
  {
    depth++;
#ifdef DEBUG
    cout << "depth: " << depth << endl;
#endif
    allleaf = true;
    vector<MultiChannelImage3DT<unsigned short int> > lastfeats = currentfeats;
    vector<vector<vector<double> > > lastRegionProbs = regionProbs;

    if ( useRegionFeature )
    {
      int a_max = ( int ) regionProbs.size();
      for ( int a = 0; a < a_max; a++ )
      {
        int b_max = ( int ) regionProbs[a].size();
        for ( int b = 0; b < b_max; b++ )
        {
          int c_max = ( int ) regionProbs[a][b].size();
          for ( int c = 0; c < c_max; c++ )
          {
            regionProbs[a][b][c] = 0.0;
          }
        }
      }
    }

#if 1
    Timer timerDepth;
    timerDepth.start();
#endif

    double weight = computeWeight ( depth, maxDepth ) - computeWeight ( depth - 1, maxDepth );

    if ( depth == 1 )
    {
      weight = computeWeight ( 1, maxDepth );
    }

//   omp_set_dynamic(0);
#pragma omp parallel for
    for ( int tree = 0; tree < nbTrees; tree++ )
    {
      const int t = ( int ) forest[tree].size();
      const int s = startnode[tree];
      startnode[tree] = t;
//#pragma omp parallel for
      for ( int i = s; i < t; i++ )
      {
        if ( !forest[tree][i].isleaf && forest[tree][i].left < 0 )
        {
          Operation *splitfeat = NULL;
          double splitval;
          double bestig = getBestSplit ( allfeats, lastfeats, labels, i, splitfeat, splitval, tree, lastRegionProbs );

          for ( int ii = 0; ii < ( int ) lastfeats.size(); ii++ )
          {
            for ( int c = 0; c < lastfeats[ii].channels(); c++ )
            {
              short unsigned int minv, maxv;
              lastfeats[ii].statistics ( minv, maxv, c );
            }
          }

          forest[tree][i].feat = splitfeat;

          forest[tree][i].decision = splitval;

          if ( splitfeat != NULL )
          {
            allleaf = 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][i].left = left;
            forest[tree][i].right = right;
            forest[tree][left].dist = vector<double> ( classes, 0.0 );
            forest[tree][right].dist = vector<double> ( classes, 0.0 );
            forest[tree][left].depth = depth;
            forest[tree][right].depth = depth;
            forest[tree][left].featcounter = 0;
            forest[tree][right].featcounter = 0;
            forest[tree][left].nodeNumber = uniquenumber;
            int leftu = uniquenumber;
            uniquenumber++;
            forest[tree][right].nodeNumber = uniquenumber;
            int rightu = uniquenumber;
            uniquenumber++;
            forest[tree][right].featcounter = 0;

#pragma omp parallel for
            for ( int iCounter = 0; iCounter < imgcounter; iCounter++ )
            {
              int xsize = currentfeats[iCounter].width();
              int ysize = currentfeats[iCounter].height();
              int zsize = currentfeats[iCounter].depth();

              for ( int x = 0; x < xsize; x++ )
              {
                for ( int y = 0; y < ysize; y++ )
                {
                  for ( int z = 0; z < zsize; z++ )
                  {
                    if ( currentfeats[iCounter].get ( x, y, z, tree ) == i )
                    {
                      Features feat;
                      feat.feats = &allfeats[iCounter];
                      feat.cfeats = &lastfeats[iCounter];
                      feat.cTree = tree;
                      feat.tree = &forest[tree];
                      feat.rProbs = &lastRegionProbs[iCounter];
                      double val = splitfeat->getVal ( feat, x, y, z );
                      if ( !isfinite ( val ) )
                      {
                        val = 0.0;
                      }

#pragma omp critical
                      if ( val < splitval )
                      {
                        currentfeats[iCounter].set ( x, y, z, left, tree );
                        if ( labelmap.find ( labels[iCounter] ( x, y, ( uint ) z ) ) != labelmap.end() )
                          forest[tree][left].dist[labelmap[labels[iCounter] ( x, y, ( uint ) z ) ]]++;
                        forest[tree][left].featcounter++;
                        if ( useCategorization && leftu < shortsize )
                          ( *globalCategorFeats[iCounter] ) [leftu]+=weight;
                      }
                      else
                      {
                        currentfeats[iCounter].set ( x, y, z, right, tree );
                        if ( labelmap.find ( labels[iCounter] ( x, y, ( uint ) z ) ) != labelmap.end() )
                          forest[tree][right].dist[labelmap[labels[iCounter] ( x, y, ( uint ) z ) ]]++;
                        forest[tree][right].featcounter++;

                        if ( useCategorization && rightu < shortsize )
                          ( *globalCategorFeats[iCounter] ) [rightu]+=weight;
                      }
                    }
                  }
                }
              }
            }

            double lcounter = 0.0, rcounter = 0.0;

            for ( uint d = 0; d < forest[tree][left].dist.size(); d++ )
            {
              if ( forbidden_classes.find ( labelmapback[d] ) != forbidden_classes.end() )
              {
                forest[tree][left].dist[d] = 0;
                forest[tree][right].dist[d] = 0;
              }
              else
              {
                forest[tree][left].dist[d] /= a[d];
                lcounter += forest[tree][left].dist[d];
                forest[tree][right].dist[d] /= a[d];
                rcounter += forest[tree][right].dist[d];
              }
            }

            if ( lcounter <= 0 || rcounter <= 0 )
            {
              cout << "lcounter : " << lcounter << " rcounter: " << rcounter << endl;
              cout << "splitval: " << splitval << " splittype: " << splitfeat->writeInfos() << endl;
              cout << "bestig: " << bestig << endl;

              for ( int iCounter = 0; iCounter < imgcounter; iCounter++ )
              {
                int xsize = currentfeats[iCounter].width();
                int ysize = currentfeats[iCounter].height();
                int zsize = currentfeats[iCounter].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 ( lastfeats[iCounter].get ( x, y, tree ) == i )
                      {
                        if ( ++counter > 30 )
                          break;

                        Features feat;

                        feat.feats = &allfeats[iCounter];
                        feat.cfeats = &lastfeats[iCounter];
                        feat.cTree = tree;
                        feat.tree = &forest[tree];
                        feat.rProbs = &lastRegionProbs[iCounter];

                        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 ( uint d = 0; d < forest[tree][left].dist.size(); d++ )
            {
              forest[tree][left].dist[d] /= lcounter;
              forest[tree][right].dist[d] /= rcounter;
            }
          }
          else
          {
            forest[tree][i].isleaf = true;
          }
        }
      }
    }


    if ( useRegionFeature )
    {
      for ( int iCounter = 0; iCounter < imgcounter; iCounter++ )
      {
        int xsize = currentfeats[iCounter].width();
        int ysize = currentfeats[iCounter].height();
        int zsize = currentfeats[iCounter].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 = currentfeats[iCounter].get ( x, y, z, tree );
                for ( uint d = 0; d < forest[tree][node].dist.size(); d++ )
                {
                  regionProbs[iCounter][ ( int ) ( allfeats[iCounter] ( x, y, z, rawChannels ) ) ][d] += forest[tree][node].dist[d];
                }
              }
            }
          }
        }
      }

      int a_max = ( int ) regionProbs.size();
      for ( int a = 0; a < a_max; a++ )
      {
        int b_max = ( int ) regionProbs[a].size();
        for ( int b = 0; b < b_max; b++ )
        {
          int c_max = ( int ) regionProbs[a][b].size();
          for ( int c = 0; c < c_max; c++ )
          {
            regionProbs[a][b][c] /= ( double ) ( rSize[a][b] );
          }
        }
      }
    }

    //compute integral images
    if ( firstiteration )
    {
      for ( int i = 0; i < imgcounter; i++ )
      {
        allfeats[i].addChannel ( ( int ) ( classes + rawChannels ) );
      }
    }

    for ( int i = 0; i < imgcounter; i++ )
    {
      computeIntegralImage ( currentfeats[i], allfeats[i], channelType.size() - classes );
    }

    if ( firstiteration )
    {
      firstiteration = false;
    }

#if 1
    timerDepth.stop();

    cout << "time for depth " << depth << ": " << timerDepth.getLastAbsolute() << endl;
#endif

    lastfeats.clear();
    lastRegionProbs.clear();
  }

  timer.stop();
  cerr << "learning finished in: " << timer.getLastAbsolute() << " seconds" << endl;
  timer.start();

  cout << "uniquenumber " << uniquenumber << endl;

  if ( useCategorization && fasthik != NULL )
  {
    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>::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 ( globalCategorFeats, ys );

  }

#ifdef DEBUG
  for ( int tree = 0; tree < nbTrees; tree++ )
  {
    int t = ( int ) forest[tree].size();

    for ( int i = 0; i < t; i++ )
    {
      printf ( "tree[%i]: left: %i, right: %i", i, forest[tree][i].left, forest[tree][i].right );

      if ( !forest[tree][i].isleaf && forest[tree][i].left != -1 )
      {
        cout <<  ", feat: " << forest[tree][i].feat->writeInfos() << " ";
        opOverview[forest[tree][i].feat->getOps() ]++;
        contextOverview[forest[tree][i].depth][ ( int ) forest[tree][i].feat->getContext() ]++;
      }

      for ( int d = 0; d < ( int ) forest[tree][i].dist.size(); d++ )
      {
        cout << " " << forest[tree][i].dist[d];
      }

      cout << endl;
    }
  }

  std::map<int, int> featTypeCounter;

  for ( int tree = 0; tree < nbTrees; tree++ )
  {
    int t = ( int ) forest[tree].size();

    for ( int i = 0; i < t; i++ )
    {
      if ( !forest[tree][i].isleaf && forest[tree][i].left != -1 )
      {
        featTypeCounter[forest[tree][i].feat->getFeatType() ] += 1;
      }
    }
  }

  cout << "evaluation of featuretypes" << endl;
  for ( map<int, int>::const_iterator it = featTypeCounter.begin(); it != featTypeCounter.end(); it++ )
  {
    cerr << it->first << ": " << it->second << endl;
  }

  for ( uint c = 0; c < ops.size(); c++ )
  {

    for ( int t = 0; t < ( int ) ops[c].size(); t++ )
    {
      cout << ops[c][t]->writeInfos() << ": " << opOverview[ops[c][t]->getOps() ] << 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 << " woContext: " << contextOverview[d][0] << " wContext: " << contextOverview[d][1] << endl;
  }
#endif

  timer.stop();
  cerr << "rest finished in: " << timer.getLastAbsolute() << " seconds" << endl;
  timer.start();
}

void SemSegContextTree::extractBasicFeatures ( NICE::MultiChannelImage3DT<double> &feats, const NICE::MultiChannelImage3DT<double> &imgData, const vector<string> &filelist, int &amountRegions )
{
  int xsize = imgData.width();
  int ysize = imgData.height();
  int zsize = imgData.depth();
  //TODO: resize image?!

  amountRegions = 0;
  feats.reInit ( xsize, ysize, zsize, imgData.channels() );
  feats.setAll ( 0 );

  for ( int z = 0; z < zsize; z++ )
  {
    NICE::MultiChannelImageT<double> feats_tmp;
    feats_tmp.reInit ( xsize, ysize, 3 );
    if ( imagetype == IMAGETYPE_RGB )
    {

      NICE::ColorImage img = imgData.getColor ( z );
      for ( int x = 0; x < xsize; x++ )
      {
        for ( int y = 0; y < ysize; y++ )
        {
          for ( int r = 0; r < 3; r++ )
          {
            feats_tmp.set ( x, y, img.getPixel ( x, y, r ), ( uint ) r );
          }
        }
      }

    }
    else
    {

      NICE::ImageT<double> img = imgData.getChannelT ( z,0 );
      for ( int x = 0; x < xsize; x++ )
      {
        for ( int y = 0; y < ysize; y++ )
        {
          feats_tmp.set ( x, y, img.getPixel ( x, y ), 0 );
        }
      }

    }

    if ( imagetype == IMAGETYPE_RGB )
      feats_tmp = ColorSpace::rgbtolab ( feats_tmp );

    for ( int x = 0; x < xsize; x++ )
    {
      for ( int y = 0; y < ysize; y++ )
      {
        if ( imagetype == IMAGETYPE_RGB )
        {
          for ( uint r = 0; r < 3; r++ )
          {
            feats.set ( x, y, z, feats_tmp.get ( x, y, r ), r );
          }
        }
        else
        {
          feats.set ( x, y, z, feats_tmp.get ( x, y, 0 ), 0 );
        }
      }
    }

    if ( useGradient )
    {
      int currentsize = feats_tmp.channels();
      feats_tmp.addChannel ( currentsize );

      for ( int c = 0; c < currentsize; c++ )
      {
        ImageT<double> tmp = feats_tmp[c];
        ImageT<double> tmp2 = feats_tmp[c+currentsize];

        NICE::FilterT<double, double, double>::gradientStrength ( tmp, tmp2 );
      }
    }

    if ( useWeijer )
    {
      if ( imagetype == IMAGETYPE_RGB )
      {
        NICE::ColorImage img = imgData.getColor ( z );
        NICE::MultiChannelImageT<double> cfeats;
        lfcw->getFeats ( img, cfeats );
        feats_tmp.addChannel ( cfeats );
      }
      else
      {
        cerr << "Can't compute weijer features of a grayscale image." << endl;
      }
    }

    if ( useGaussian )
    {
      vector<string> list;
      StringTools::split ( filelist[z], '/', list );
      string gaussPath = StringTools::trim ( filelist[z], list.back() ) + "gaussmap/" + list.back();
      NICE::Image gauss ( gaussPath );
      feats_tmp.addChannel ( gauss );
      //cout << "Added file " << gaussPath << " to feature stack " << endl;
    }
    
    // read the geometric cues produced by Hoiem et al.
    if ( useHoiemFeatures )
    {
      // 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 );

      // Now we have to do some fancy regular expressions :)
      // Original image filename: basel_000083.jpg
      // hoiem result: basel_000083_c_sky.png

      // Fancy class of Ferid which supports string handling especially for filenames
      FileName fn ( filelist[z] );
      fn.removeExtension();
      FileName fnBase = fn.extractFileName();

      // counter for the channel index, starts with the current size of the destination multi-channel image
      int currentChannel = feats_tmp.channels();

      // add a channel for each feature in advance
      feats_tmp.addChannel ( hoiemClasses.size() );

      // loop through all geometric categories and add the images
      for ( vector<string>::const_iterator i = hoiemClasses.begin(); i != hoiemClasses.end(); i++, currentChannel++ )
      {
        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() );
          // check whether the image size is consistent
          if ( confidenceImage.width() != feats_tmp.width() || confidenceImage.height() != feats_tmp.height() )
          {
            fthrow ( Exception, "The size of the geometric confidence image does not match with the original image size: " << fnConfidenceImage.str() );
          }
          ImageT<double> dst = feats_tmp[currentChannel];

          // copy standard image to double image
          for ( uint y = 0 ; y < ( uint ) confidenceImage.height(); y++ )
            for ( uint x = 0 ; x < ( uint ) confidenceImage.width(); x++ )
              feats_tmp ( x, y, currentChannel ) = ( double ) confidenceImage ( x, y );
        }
      }
    }

    uint oldChannels = feats.channels();
    if ( feats.channels() < feats_tmp.channels() )
      feats.addChannel ( feats_tmp.channels()-feats.channels() );

    for ( int x = 0; x < xsize; x++ )
    {
      for ( int y = 0; y < ysize; y++ )
      {
        for ( uint r = oldChannels; r < ( uint ) feats_tmp.channels(); r++ )
        {
          feats.set ( x, y, z, feats_tmp.get ( x, y, r ), r );
        }
      }
    }
  }

  if ( useRegionFeature )
  {
    //using segmentation
    MultiChannelImageT<int> regions;
    regions.reInit( xsize, ysize, zsize );
    vector<int> chanSelect;
    for ( int i=0; i<3; i++ )
      chanSelect.push_back ( i );

    amountRegions = segmentation->segRegions ( imgData, regions, chanSelect );

    int cchannel = feats.channels();
    feats.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++ )
        {
          feats.set ( x, y, z, regions ( x, y, ( uint ) z ), cchannel );
        }
      }
    }
  }
  
  if ( useVariance )
  {
    int cchannel = feats.channels();
    if (imagetype = IMAGETYPE_RGB)
    {
      feats.addChannel( 3 );
      for (int c = 0; c < 3; c++)
      {
        feats.calcVariance( c, cchannel+c );
      }
    }
    else
    {
      feats.addChannel( 1 );
      feats.calcVariance( 0, cchannel );
    }
  }

}

void SemSegContextTree::semanticseg ( 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();

  firstiteration = true;

  int classes = labelmapback.size();

  int numClasses = classNames->numClasses();

  fprintf ( stderr, "ContextTree classification !\n" );

  probabilities.reInit ( xsize, ysize, zsize, numClasses );
  probabilities.setAll ( 0 );

  SparseVector *globalCategorFeat = new SparseVector();

  MultiChannelImage3DT<double> feats;

  // Basic Features
  int amountRegions;
  extractBasicFeatures ( feats, imgData, filelist, amountRegions ); //read image and do some simple transformations

  vector<int> rSize;
  if ( useRegionFeature )
  {
    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[feats ( x, y, z, rawChannels ) ]++;
        }
      }
    }
  }

  bool allleaf = false;

  MultiChannelImage3DT<unsigned short int> currentfeats ( xsize, ysize, zsize, nbTrees );

  currentfeats.setAll ( 0 );

  depth = 0;

  vector<vector<double> > regionProbs;
  if ( useRegionFeature )
  {
    regionProbs = vector<vector<double> > ( amountRegions, vector<double> ( classes, 0.0 ) );
  }

  for ( int d = 0; d < maxDepth && !allleaf; d++ )
  {
    depth++;
    vector<vector<double> > lastRegionProbs = regionProbs;
    if ( useRegionFeature )
    {
      int b_max = ( int ) regionProbs.size();
      for ( int b = 0; b < b_max; b++ )
      {
        int c_max = ( int ) regionProbs[b].size();
        for ( int c = 0; c < c_max; c++ )
        {
          regionProbs[b][c] = 0.0;
        }
      }
    }

    double weight = computeWeight ( depth, maxDepth ) - computeWeight ( depth - 1, maxDepth );

    if ( depth == 1 )
    {
      weight = computeWeight ( 1, maxDepth );
    }

    allleaf = true;

    MultiChannelImage3DT<unsigned short int> lastfeats = currentfeats;

    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 t = currentfeats.get ( x, y, z, tree );

            if ( forest[tree][t].left > 0 )
            {
              allleaf = false;
              Features feat;
              feat.feats = &feats;
              feat.cfeats = &lastfeats;
              feat.cTree = tree;
              feat.tree = &forest[tree];
              feat.rProbs = &lastRegionProbs;

              double val = forest[tree][t].feat->getVal ( feat, x, y, z );
              if ( !isfinite ( val ) )
              {
                val = 0.0;
              }

              if ( val < forest[tree][t].decision )
              {
                currentfeats.set ( x, y, z, forest[tree][t].left, tree );
#pragma omp critical
                {
                  if ( fasthik != NULL && useCategorization && forest[tree][forest[tree][t].left].nodeNumber < uniquenumber )
                    ( *globalCategorFeat ) [forest[tree][forest[tree][t].left].nodeNumber] += weight;
                }
              }
              else
              {
                currentfeats.set ( x, y, z, forest[tree][t].right, tree );
#pragma omp critical
                {
                  if ( fasthik != NULL && useCategorization && forest[tree][forest[tree][t].right].nodeNumber < uniquenumber )
                    ( *globalCategorFeat ) [forest[tree][forest[tree][t].right].nodeNumber] += weight;
                }
              }
            }
          }
        }
      }
    }

    if ( useRegionFeature )
    {
      int xsize = currentfeats.width();
      int ysize = currentfeats.height();
      int zsize = currentfeats.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 = currentfeats.get ( x, y, z, tree );
              for ( uint d = 0; d < forest[tree][node].dist.size(); d++ )
              {
                regionProbs[ ( int ) ( feats ( x, y, z, rawChannels ) ) ][d] += forest[tree][node].dist[d];
              }
            }
          }
        }
      }

      for ( int b = 0; b < ( int ) regionProbs.size(); b++ )
      {
        for ( int c = 0; c < ( int ) regionProbs[b].size(); c++ )
        {
          regionProbs[b][c] /= ( double ) ( rSize[b] );
        }
      }
    }

    if ( depth < maxDepth )
    {
      //compute integral images
      if ( firstiteration )
      {
        feats.addChannel ( classes + rawChannels );
      }
      computeIntegralImage ( currentfeats, feats, channelType.size() - classes );
      if ( firstiteration )
      {
        firstiteration = false;
      }
    }
  }

  int allClasses = ( int ) probabilities.channels();
  vector<int> useclass ( allClasses, 1 );

  vector<int> classesInImg;

  if ( useCategorization )
  {
    if ( cndir != "" )
    {
      for ( int z = 0; z < zsize; z++ )
      {
        std::vector< std::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 < allClasses );
          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 );
    }
  }

  if ( pixelWiseLabeling )
  {
    //finales labeln:
    //long int offset = 0;

    for ( int x = 0; x < xsize; x++ )
    {
      for ( int y = 0; y < ysize; y++ )
      {
        for ( int z = 0; z < zsize; z++ )
        {
          double maxvalue = - numeric_limits<double>::max(); //TODO: das kann auch nur pro knoten gemacht werden, nicht pro pixel
          int maxindex = 0;

          for ( uint c = 0; c < classesInImg.size(); c++ )
          {
            int i = classesInImg[c];
            int currentclass = labelmapback[i];
            if ( useclass[currentclass] )
            {
              probabilities ( x, y, z, currentclass ) = getMeanProb ( x, y, z, i, currentfeats );

              if ( probabilities ( x, y, z, currentclass ) > maxvalue )
              {
                maxvalue = probabilities ( x, y, z, currentclass );
                maxindex = currentclass;
              }
            }
          }
          segresult.set ( x, y, maxindex, ( uint ) z );
          if ( maxvalue > 1 )
            cout << "maxvalue: " << maxvalue << endl;
        }
      }
    }
#undef VISUALIZE
#ifdef VISUALIZE
    for ( int z = 0; z < zsize; z++ )
    {
      for ( int j = 0 ; j < ( int ) probabilities.numChannels; j++ )
      {
        //cout << "class: " << j << endl;//" " << cn.text (j) << endl;

        NICE::Matrix tmp ( probabilities.height(), probabilities.width() );
        double maxval = -numeric_limits<double>::max();
        double minval = numeric_limits<double>::max();


        for ( int y = 0; y < probabilities.height(); y++ )
          for ( int x = 0; x < probabilities.width(); x++ )
          {
            double val = probabilities ( x, y, z, j );
            tmp ( y, x ) = val;
            maxval = std::max ( val, maxval );
            minval = std::min ( val, minval );
          }
        tmp ( 0, 0 ) = 1.0;
        tmp ( 0, 1 ) = 0.0;

        NICE::ColorImage imgrgb ( probabilities.width(), probabilities.height() );
        ICETools::convertToRGB ( tmp, imgrgb );

        cout << "maxval = " << maxval << " minval: " << minval << " for class " << j << endl; //cn.text (j) << endl;

        std::string s;
        std::stringstream out;
        out << "tmpprebmap" << j << ".ppm";
        s = out.str();
        imgrgb.write ( s );
        //showImage(imgrgb, "Ergebnis");
        //getchar();
      }
    }
    cout << "fertsch" << endl;
    getchar();
    cout << "weiter gehtsch" << endl;
#endif
  }
  else
  {
    //using segmentation
    NICE::MultiChannelImageT<int> regions;
    int xsize = feats.width();
    int ysize = feats.height();
    int zsize = feats.depth();
    regions.reInit ( xsize, ysize, zsize );

    if ( useRegionFeature )
    {
      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, feats ( x, y, z, rchannel ), ( uint ) z );
          }
        }
      }
    }
    else
    {
      amountRegions = 0;
      vector<int> chanSelect;
      for ( int i=0; i<3; i++ )
        chanSelect.push_back ( i );
      amountRegions = segmentation->segRegions ( imgData, regions, chanSelect );

#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 cregion = regions ( x, y, ( uint ) z );
          for ( uint c = 0; c < classesInImg.size(); c++ )
          {
            int d = classesInImg[c];
            regionProbs[cregion][d] += getMeanProb ( x, y, z, d, currentfeats );
          }
        }
      }
    }

    for ( int r = 0; r < amountRegions; r++ )
    {
      double maxval = regionProbs[r][classesInImg[0]];
      bestlabels[r] = classesInImg[0];

      for ( int d = 1; d < classes; d++ )
      {
        if ( maxval < regionProbs[r][d] )
        {
          maxval = regionProbs[r][d];
          bestlabels[r] = d;
        }
      }

      bestlabels[r] = labelmapback[bestlabels[r]];
    }
    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 );
        }
      }
    }

//#define WRITEREGIONS
#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
  }

  cout << "segmentation finished" << endl;
}

void SemSegContextTree::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 << integralMap.size() << endl;
  for ( int i = 0; i < ( int ) integralMap.size(); i++ )
  {
    os << integralMap[i].first << " " << integralMap[i].second << endl;
  }

  os << rawChannels << endl;

  os << uniquenumber << endl;
}

void SemSegContextTree::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 )
            {
              for ( uint c = 0; c < ops[o].size(); c++ )
              {
                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;
    channelType.push_back ( tmp );
  }

  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;
}