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

  // 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 *= 3;
    }
    
    if(ft == 1)
    {
      if(depth < 8)
      {
        ft = 0;
      }
    }

		// random value range between (-window-size/2) and (window-size/2)
		// TODO z-values either positive or negative
		int tmp_z =  feats[0].depth()/6;
    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);
    }
    
    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;

  // 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();
		NICE::MultiChannelImageT<uchar> pixelLabels;
		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;
		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);
				}
			}
		}
	}

}

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