//
// This file is part of 3DCellAtlas. 
// Copyright (C) 2015 George W. Bassel and collaborators.
//
// If you use 3DCellAtlas in your work, please cite:
//   http://dx.doi.org/10.1105/tpc.15.00175
//
// 3DCellAtlas is an AddOn for MorphoGraphX - http://www.MorphoGraphX.org
// Copyright (C) 2012-2015 Richard S. Smith and collaborators.
//
// 3DCellAtlas and MorphoGraphX are free software, and are licensed under under the terms of the 
// GNU General (GPL) Public License version 2.0, http://www.gnu.org/licenses.
// 
#include "RootCellProcessing.hpp"

#include "Triangulate.hpp"

// gsl for splines
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>

#include <gsl/gsl_bspline.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics.h>
namespace mgx {

  void writeAttrMaps(RootCellProcessing rcp, CellAtlasAttr *cellAtlasAttr, CellAtlasConfigAttr *cellAtlasConfigAttr)
  {

    (*cellAtlasConfigAttr)[0].wallArea = rcp.rootData.wallArea;
    (*cellAtlasConfigAttr)[0].diffBezInterp = rcp.rootData.diffBezInterp;
    (*cellAtlasConfigAttr)[0].bez = rcp.rootData.bez;
    //(*cellAtlasConfigAttr)[0].bezLength = rcp.rootData.bezLength;
    (*cellAtlasConfigAttr)[0].analyzed = true;

    //cout << "bla " << rcp.rootData.bez.size() << "/" << rcp.rootData.diffBezInterp.size() << endl;

    for(int i = 0; i<rcp.rootData.numCells; i++){
      int currentLabel = rcp.rootData.uniqueLabels[i];
      (*cellAtlasAttr)[i].cellLabel = currentLabel;
      (*cellAtlasAttr)[i].dirRad = rcp.rootData.dirRad[currentLabel];
      (*cellAtlasAttr)[i].dirLong = rcp.rootData.dirLong[currentLabel];
      (*cellAtlasAttr)[i].dirCirc = rcp.rootData.dirCirc[currentLabel];

      (*cellAtlasAttr)[i].sizeRad = rcp.rootData.lengthRad[currentLabel];
      (*cellAtlasAttr)[i].sizeLong = rcp.rootData.lengthLong[currentLabel];
      (*cellAtlasAttr)[i].sizeCirc = rcp.rootData.lengthCirc[currentLabel];

      (*cellAtlasAttr)[i].coordLong = rcp.rootData.arclengths[currentLabel];
      (*cellAtlasAttr)[i].coordLongUM = rcp.rootData.arclengthsUM[currentLabel];

      (*cellAtlasAttr)[i].coordLongMinUM = rcp.rootData.arclengthsMinUM[currentLabel];
      (*cellAtlasAttr)[i].coordLongMaxUM = rcp.rootData.arclengthsMaxUM[currentLabel];

      (*cellAtlasAttr)[i].nrNeighbors = rcp.rootData.nrNeighbors[currentLabel];

      (*cellAtlasAttr)[i].coordCirc = rcp.rootData.azimuthal[currentLabel];
      (*cellAtlasAttr)[i].coordRad = rcp.rootData.scaledRadialDis[currentLabel];
      (*cellAtlasAttr)[i].coordRadAbs = rcp.rootData.radialDis[currentLabel];

      (*cellAtlasAttr)[i].centroid = rcp.rootData.cellCentroids[currentLabel];
      (*cellAtlasAttr)[i].volume = rcp.rootData.cellVolumes[currentLabel];
      (*cellAtlasAttr)[i].signalCentroid = rcp.rootData.signalCentroids[currentLabel];

     // (*cellAtlasAttr)[i].signalCentroid = rcp.rootData.signalCentroids[currentLabel];
     // (*cellAtlasAttr)[i].signalCentroid = rcp.rootData.signalCentroids[currentLabel];

      (*cellAtlasAttr)[i].nearestSurfacePoint = rcp.rootData.nearestSurfacePoint[currentLabel];
      (*cellAtlasAttr)[i].wallArea = rcp.rootData.cellWallArea[currentLabel];
      (*cellAtlasAttr)[i].outsideWallArea = rcp.rootData.outsideWallArea[currentLabel];
      (*cellAtlasAttr)[i].outsideWallAreaPercent = rcp.rootData.outsideWallArea[currentLabel]/rcp.rootData.cellWallArea[currentLabel];
      (*cellAtlasAttr)[i].associatedCorticalCell = rcp.rootData.associatedCorticalCell[currentLabel];

      (*cellAtlasAttr)[i].badCell = rcp.rootData.cellBad[currentLabel];
    }

    for(int i = 0; i < rcp.rootData.numCells; i++){
      (*cellAtlasConfigAttr)[0].labelIdxMap[(*cellAtlasAttr)[i].cellLabel] = i;
    }
    

  }

// convert the attr map information into a single map (label->data)
std::map<int, double> RootCellProcessing::getMapFromAttributes(QString choice)
{
  std::map<int, double> result;

  //int numCells = (*data).size();

  for(int i = 0; i<numCells; i++){
    int label = labelMap[i];//(*data)[i].cellLabel;
    double value = 0;
    if(choice == "Longitudinal (Arclengths)"){
      value = (*data)[i].coordLong;
    } else if(choice == "Radial"){
      value = (*data)[i].coordRad;
    } else if(choice == "Radial (Absolute)"){
      value = (*data)[i].coordRadAbs;
    } else if(choice == "Longitudinal Cell Length"){
      value = (*data)[i].sizeLong;
    } else if(choice == "Radial Cell Length"){
      value = (*data)[i].sizeRad;
    } else if(choice == "Circumferential Cell Length"){
      value = (*data)[i].sizeCirc;
    } else if(choice == "Volume"){
      value = (*data)[i].volume;
    } else if(choice == "Circumferential Angle"){
      value = (*data)[i].coordCirc;
    } else if(choice == "Cell Wall Area"){
      value = (*data)[i].wallArea;
    } else if(choice == "Associated Cortical Cell"){
      value = (*data)[i].associatedCorticalCell;
    } else if(choice == "Cell Volume"){
      value = (*data)[i].volume;
    } else if(choice == "Cell Wall Area"){
      value = (*data)[i].wallArea;
    }

    result[label] = value;
  }

  return result;

}



void RootCellProcessing::initHeatMapDataStructure(heatMapDataStructure& body, double sig, int& chosenX, int& chosenY)
{
  body.gridSize = 160;
  body.sigma = sig;

  QString option = "Circumferential Angle";

  if(chosenX == 2){
    option = "Circumferential Angle";//rootData.azimuthal;
  } else if(chosenX == 3) {
    option = "Longitudinal (Arclengths)";//rootData.arclengths;
  } else if(chosenX == 4) {
    option = "Radial Cell Length";//rootData.lengthRad;
  } else if(chosenX == 5) {
    option = "Circumferential Cell Length";//rootData.lengthCirc;
  } else if(chosenX == 6) {
    option = "Longitudinal Cell Length";//rootData.lengthLong;
  } else if(chosenX == 7) {
    option = "Cell Volume";//rootData.lengthLong;
  } else if(chosenX == 8) {
    option = "Cell Wall Area";//rootData.lengthLong;
  } else {
    option = "Radial";//rootData.scaledRadialDis;
    chosenX = 1;
  }

  body.x = getMapFromAttributes(option);

  if(chosenY == 1){
    option = "Radial";//rootData.scaledRadialDis;
  } else if(chosenY == 2) {
    option = "Circumferential Angle";//rootData.azimuthal;
  } else if(chosenY == 3) {
    option = "Longitudinal (Arclengths)";//rootData.arclengths;
  } else if(chosenY == 4) {
    option = "Radial Cell Length";//rootData.lengthRad;
  } else if(chosenY == 6) {
    option = "Longitudinal Cell Length";//rootData.lengthLong;
  } else if(chosenY == 7) {
    option = "Cell Volume";//rootData.lengthLong;
  } else if(chosenY == 8) {
    option = "Cell Wall Area";//rootData.lengthLong;
  } else {
    option = "Circumferential Cell Length";//rootData.lengthCirc;
    chosenY = 5;
  }

  body.y = getMapFromAttributes(option);
}

// finds min, max, avg and std of the cells, needed for scaling the heatmap size
void RootCellProcessing::findHeatStatistics(heatMapDataStructure& body)
{

  double minX = 1E20, minY = 1E20, maxX = -1E20, maxY = -1E20;
  double avgX = 0, avgY = 0, stdX = 0, stdY = 0;

  //int numCells = (*data).size();

  for(int i = 0; i<numCells; i++){
    int currentLabel = labelMap[i];//(*data)[i].cellLabel;
    if(!body.preselect[currentLabel] and body.cells.find(currentLabel) != body.cells.end()){
      // NAN check
      if(body.x[currentLabel]!=body.x[currentLabel] || body.y[currentLabel]!=body.y[currentLabel]){
        body.x[currentLabel] = 0;
        body.y[currentLabel] = 0;
        cout << " * * * Warning: Cell with label " << currentLabel << " was detected with invalid cell sizes and set to size zero" << endl; 
      }

      avgX += body.x[currentLabel];
      avgY += body.y[currentLabel];
      if(body.cells.find(currentLabel) != body.cells.end()){
      //cout << "lab  " << currentLabel << "  " << x[currentLabel] << "  Y " << y[currentLabel] << endl;
        if(minX > body.x[currentLabel]){
          minX = body.x[currentLabel];}
        if(minY > body.y[currentLabel]){
          minY = body.y[currentLabel];}
        if(maxX < body.x[currentLabel]){
          maxX = body.x[currentLabel];}
        if(maxY < body.y[currentLabel]){
//cout << "newMaxY  " << currentLabel << "  " << body.x[currentLabel] << "  Y " << body.y[currentLabel] << endl;
          maxY = body.y[currentLabel];}
      }
    }
  }

  avgX /= body.cellCounter;
  avgY /= body.cellCounter;

  for(int i = 0; i<numCells; i++){
    int currentLabel = labelMap[i];//(*data)[i].cellLabel;
    if(!body.preselect[currentLabel] and body.cells.find(currentLabel) != body.cells.end()){
      stdX += std::abs(avgX - body.x[currentLabel]);
      stdY += std::abs(avgY - body.y[currentLabel]);
    }
  }

  stdX /= body.cellCounter;
  stdY /= body.cellCounter;

  //cout << "minmax  " << minX << "  " << maxX << " / " << minY << " " << maxY << endl;

  body.heatMapMinMax[0] = minX;
  body.heatMapMinMax[1] = maxX;
  body.heatMapMinMax[2] = minY;
  body.heatMapMinMax[3] = maxY;

  body.heatMapStatistics[0] = avgX;
  body.heatMapStatistics[1] = avgY;
  body.heatMapStatistics[2] = stdX;
  body.heatMapStatistics[3] = stdY;

}

// generates the heatmap (=cellMorphLandscape) of the cell set (=body)
cellMorphLandscape RootCellProcessing::generateCellMorphologyLandscape(heatMapDataStructure& body)
{

  cellMorphLandscape xyLandscape;

  findHeatStatistics(body);

  //int numCells = (*data).size();

  double minX = body.heatMapMinMax[0];
  double minY = body.heatMapMinMax[2];
  double maxX = body.heatMapMinMax[1];
  double maxY = body.heatMapMinMax[3];

  double avgX = body.heatMapStatistics[0];
  double avgY = body.heatMapStatistics[1];
  double stdX = body.heatMapStatistics[2];
  double stdY = body.heatMapStatistics[3];

  //cout << " mm     " << minX << "  " << maxX << "  " << minY << "  " << maxY << endl;
  //cout << " stat   " << avgX << "  " << stdX << "  " << avgY << "  " << stdY << endl;

  if(outlierDetection){
    double threshold = 20; // threshold for outlier detection

    // take care of outliers
    if(avgX + threshold*stdX < maxX or avgX - threshold*stdX > minX or avgY + threshold*stdY < maxY or avgY - threshold*stdY > minY){
      // remove outlier cells and rerun statistics
      for(int i = 0; i<numCells; i++){
        int currentLabel = labelMap[i];//(*data)[i].cellLabel;
        if(!body.preselect[currentLabel] and body.cells.find(currentLabel) != body.cells.end()){
          if(avgX + threshold*stdX < body.x[currentLabel] or avgX - threshold*stdX > body.x[currentLabel] or 
             avgY + threshold*stdY < body.y[currentLabel] or avgY - threshold*stdY > body.y[currentLabel]){
            // this might cause problems with the preselect feature!
            body.x[currentLabel] = 0;
            body.y[currentLabel] = 0;
            cout << " * * * Warning: Cell with label " << currentLabel << " was detected as an outlier with abnormal cell sizes and set to size zero" << endl; 
          }
        }
      }
      findHeatStatistics(body);
      minX = body.heatMapMinMax[0];
      minY = body.heatMapMinMax[2];
      maxX = body.heatMapMinMax[1];
      maxY = body.heatMapMinMax[3];
    }
  }
  // step size of the heatmap
  double stepX = (double)(maxX-minX)/(double)body.gridSize;
  double stepY = (double)(maxY-minY)/(double)body.gridSize;

  for(int i = 0; i<body.gridSize; i++){
    for(int j = 0; j<body.gridSize; j++){
      xyLandscape[i][j] = 0;
      body.points[i][j] = 0;
    }
  }

  for(int i = 0; i<numCells; i++){
    // calc center of gaussian
    int currentLabel = labelMap[i];//(*data)[i].cellLabel;
    if(!body.preselect[currentLabel] and body.cells.find(currentLabel) != body.cells.end()){
      int binX = interpolateArrayIndex(body.x[currentLabel],minX,maxX,body.gridSize-1);
      int binY = interpolateArrayIndex(body.y[currentLabel],minY,maxY,body.gridSize-1);

      body.points[binX][binY]++;
      body.pointsPos[currentLabel] = Point2d(binX,binY);

      // cut off radius around the center
      int radius = std::max(10.0,body.sigma*2);
      for(int i = -radius; i<=radius; i++){
        for(int j = -radius; j<=radius; j++){
          // calc gaussian
          if(binX+i>=0 and binX+i<body.gridSize and binY+j>=0 and binY+j<body.gridSize){
            double xValue = interpolateArrayIndex(body.x[currentLabel]+i*stepX,minX,maxX,body.gridSize-1);
            double yValue = interpolateArrayIndex(body.y[currentLabel]+j*stepY,minY,maxY,body.gridSize-1);
            xyLandscape[binX+i][binY+j] += gauss2D(xValue, yValue, binX, binY, body.sigma, body.sigma);
          }
        }
      }
    }
  }
  //cout << " generated     " << endl;
  return xyLandscape;
}

// find all local maxima in a cellMorphLandscape (=heatmap) and save their positions
std::map<int,Point2d> RootCellProcessing::findMaxima(cellMorphLandscape& landsc)
{
  std::map<int,Point2d> maxima;

  int counter = 0;
  int size = landsc.size();

  // check every point of the landscape
  for(int i=0; i<size; i++){
    for(int j=0; j<size; j++){
      bool isMax = true;
      for(int in=-1; in<=1; in++){
        for(int jn=-1; jn<=1; jn++){
          if(i+in >= 0 and i+in < size and j+jn >= 0 and j+jn < size and (jn!=0 or in!=0))
            if(landsc[i][j] <= landsc[i+in][j+jn]) isMax = false;
        }
      }
      if(isMax){
        Point2d maxPos(i,j);
        maxima[counter] = maxPos;
        counter++;
      }
    }
  }

  return maxima;
  }

// label all cells according to a label-to-maximum and a maximum-to-referenceMaximum map
std::map<int,double> RootCellProcessing::labelCells(std::map<int,double>& labelToMax, std::map<int,double>& maxToRef, heatMapDataStructure& body)
{
  std::set<int> rootArea = body.cells;
  IntIntMap labels = body.maximaHeatMapSelectedLabel;


  std::map<int,double> cellLabel;
  // generate data for the graphs

  //int numCells = (*data).size();

  for(int i = 0; i<numCells; i++){
    int currentLabel = labelMap[i];//(*data)[i].cellLabel;
    if(!body.preselect[currentLabel] and body.cells.find(currentLabel) != body.cells.end()){
      cellLabel[currentLabel] = labels[maxToRef[labelToMax[currentLabel]]];
      //cout << "label cells " << currentLabel << "/" << cellLabel[currentLabel] << endl;
    }
  }
  return cellLabel;
  }

// reverse the arclengths if the root has the wrong orientation
void RootCellProcessing::correctDirections()
  {

  // check azimut of first cell

  //int numCells = (*data).size();

  int idxFirst = (*config)[0].labelIdxMap[(*config)[0].labelFirstCell];

  if((*data)[idxFirst].coordLong < 0.5) {
  // if low -> right direction
  // do nothing
  } else {
    // if high -> wrong direction
    // reverse arclengths
    double bezLength = 0;

    for(int j=1; j<(*config)[0].bez.size(); j++){
      bezLength += norm((*config)[0].bez[j]-(*config)[0].bez[j-1]);
    }

    for(int i = 0; i<numCells; i++){
      (*data)[i].coordLong = 1.0 - (*data)[i].coordLong;
      (*data)[i].coordLongUM = bezLength - (*data)[i].coordLongUM;
    }

    int save = (*config)[0].labelLastRootCapCell;
    (*config)[0].labelLastRootCapCell = (*config)[0].labelFirstCell;
    (*config)[0].labelFirstCell = save;

    assignRootRegions();
  }

  }

// assign the root regions according the selected cells and arclengths
void RootCellProcessing::assignRootRegions()
  {

  //int numCells = (*data).size();

  mainBody.cellCounter = 0;
  columella.cellCounter = 0;
  radicle.cellCounter = 0;

  std::set<int> emptySet;

  mainBody.cells = emptySet;
  columella.cells = emptySet;
  radicle.cells = emptySet;

  if(!(*config)[0].twoParts){
    // no root cap -> all cells are in the main body
    for(int i=0; i<numCells; i++){
      int currentLabel = labelMap[i];//(*data)[i].cellLabel;
      mainBody.cells.insert(currentLabel);
      mainBody.cellCounter++;
    }
  } else {
    // with root cap (two parts) -> 
    int idxFirst = (*config)[0].labelIdxMap[(*config)[0].labelFirstCell];
    int idxLastRootCap = (*config)[0].labelIdxMap[(*config)[0].labelLastRootCapCell];
    double arcFirst = (*data)[idxFirst].coordLong;
    double arcLastRootCap = (*data)[idxLastRootCap].coordLong;

    for(int i=0; i<numCells; i++){
      (*data)[i].isPart1 = false;
      (*data)[i].isPart2 = false;
      (*data)[i].isColumella = false;
      double currentLong = (*data)[i].coordLong;
      int currentLabel = labelMap[i];//(*data)[i].cellLabel;
      if(currentLong <= arcFirst){
        columella.cells.insert(currentLabel);
        (*data)[i].isColumella = true;
        columella.cellCounter++;
      } else if(currentLong <= arcLastRootCap){
        radicle.cells.insert(currentLabel);
        (*data)[i].isPart2 = true;
        radicle.cellCounter++;
      } else {
        mainBody.cells.insert(currentLabel);
        (*data)[i].isPart1 = true;
        mainBody.cellCounter++;
      }
    }
  }

}

// relate each cell in the root to a maxima in the heatmap by following the steepest gradient upwards until a maximum is reached
std::map<int,double> RootCellProcessing::relateCellstoMaxima(heatMapDataStructure& body, cellMorphLandscape& landsc, std::map<int,double>& x, 
  std::map<int,double>& y, std::map<int,Point2d> maxima)
{
  std::map<int,double> nearestMaximum;
  std::set<int> rootArea = body.cells;
  int numMax = maxima.size();
  double minX = body.heatMapMinMax[0];
  double maxX = body.heatMapMinMax[1];
  double minY = body.heatMapMinMax[2];
  double maxY = body.heatMapMinMax[3];

  //int numCells = (*data).size();

  for(int i = 0; i<numCells; i++){
    int currentLabel = labelMap[i];//(*data)[i].cellLabel;//rootData.uniqueLabels[i];
    if(rootArea.find(currentLabel) != rootArea.end()){
      Point2d currentPos;
      double minDis = 1E20;
      int minNum = 0;
      currentPos[0] = interpolateArrayIndex(x[currentLabel],minX,maxX,body.gridSize-1);
      currentPos[1] = interpolateArrayIndex(y[currentLabel],minY,maxY,body.gridSize-1);
      for(int j = 0; j<numMax; j++){ // find nearest maximum
        if(minDis > norm(currentPos - maxima[j])){
          minDis = norm(currentPos - maxima[j]);
          minNum = j;
        }
      }

      int stopCounter = 0;

      // if nearest maximum too far, then follow gradient
      while(minDis > body.sigma and stopCounter<100){
        stopCounter++;
        double maxValue = -1E20;
        int changeX = 0;
        int changeY = 0;
        // find highest value in neighborhood and go there
        for(int i = -1; i<=1; i++){
          for(int j = -1; j<=1; j++){
            if(currentPos[0]+i >= 0 and currentPos[0]+i < body.gridSize and currentPos[1]+j >=0 and currentPos[1]+j < body.gridSize){
              if(maxValue < landsc[currentPos[0]+i][currentPos[1]+j]){
                maxValue = landsc[currentPos[0]+i][currentPos[1]+j];
                changeX = i;
                changeY = j;
              }
            }
          }
        }
        // find nearest max again
        minDis = 1E20;
        minNum = 0;
        currentPos[0] += changeX;
        currentPos[1] += changeY;
        for(int j = 0; j<numMax; j++){ // find nearest maximum
          if(minDis > norm(currentPos - maxima[j])){
            minDis = norm(currentPos - maxima[j]);
            minNum = j;
          }
        }
      }
      nearestMaximum[currentLabel] = minNum;
    }
  }


  return nearestMaximum;
}

// calculates the mininmum distance to a maximum from a given maximum
double RootCellProcessing::minDisToMax(std::map<int,Point2d>& maxima, std::map<int,double>& usedMax, int currentIdx)
{
  int numMax = maxima.size();
  double minDistance = 1E20;
  for(int j = 0; j<numMax; j++){
    if(usedMax[j] == 1){
      double distance = norm(maxima[j] - maxima[currentIdx]);
      if(distance < minDistance){
        minDistance = distance;
      }
    }
  }
  return minDistance;
}

// function to assign cluster automatically by selecting the highest maxima
void RootCellProcessing::findAutomaticCluster(heatMapDataStructure& body, int nrOfClusters)
{
  cellMorphLandscape landsc = body.heatMap;
  std::map<int,Point2d> maxima = body.maximaHeatMapAll;

  std::map<int,Point2d> automaticCluster;
  std::map<int,double> valueOfMax;
  std::map<int,double> usedMax;
  int numMax = maxima.size();
  IntIntMap labels;

  for(int j = 0; j<numMax; j++){
    Point2d currentPoint = maxima[j];
    valueOfMax[j] = landsc[currentPoint[0]][currentPoint[1]];
    usedMax[j] = 0;
  }


  for(int i = 0; i<nrOfClusters; i++){
    int maxValue = 0;
    int maxPos = -1;
    for(int j = 0; j<numMax; j++){
      if(usedMax[j] == 0){
        if(maxValue < valueOfMax[j] && minDisToMax(maxima,usedMax,j) > 10.0){
          maxValue = valueOfMax[j];
          maxPos = j;
        }
      }
    }
    usedMax[maxPos] = 1;
    automaticCluster[i] = maxima[maxPos];
    labels[i] = i;
  }

  body.maximaHeatMapSelected = automaticCluster;
  body.maximaHeatMapSelectedLabel = labels;

}

void RootCellProcessing::setAutoCluster(int nrCluster, bool ismainBody)
{
  if(ismainBody)
    findAutomaticCluster(mainBody, nrCluster);
  else
    findAutomaticCluster(radicle, nrCluster);

}

// deletes a cluster from the current selected ones
void RootCellProcessing::delIdx(heatMapDataStructure& body, int idx)
{

  std::map<int,Point2d> max = body.maximaHeatMapSelected;
  IntIntMap labels = body.maximaHeatMapSelectedLabel;

  int sizeMax = max.size();
  std::map<int,Point2d> newMax;
  IntIntMap newLabels;
  int correct = 0;

  for(int i = 0; i<sizeMax; i++){
    if(i!=idx){
      newMax[i-correct] = max[i];
      newLabels[i-correct] = labels[i];
    } else {
      correct++;
    }
  }
  body.maximaHeatMapSelected = newMax;
  body.maximaHeatMapSelectedLabel = newLabels;
}


// big procedure to find the correct labels, calls earlier functions
std::map<int,double> RootCellProcessing::clusterCells(heatMapDataStructure& body)
{

  cellMorphLandscape landsc = body.heatMap;
  std::map<int,double> x = body.x;
  std::map<int,double> y = body.y;
  std::map<int,Point2d> maxima = body.maximaHeatMapAll;
  std::map<int,Point2d> customMax = body.maximaHeatMapSelected;

  // relate each cell to the nearest maximum
  std::map<int,double> nearestMaximum = relateCellstoMaxima(body, landsc, x, y, maxima);

  int numMax = maxima.size();
  int nrOfCustomMax = customMax.size();

  std::map<int,double> maximaClustMap;

  for(int i = 0; i<numMax; i++){
    // calc distance to reference maximum and find nearest ref max
    double minDis = 1E20;
    int idxRef = 0;
    for(int j = 0; j<nrOfCustomMax; j++){
      double distance = norm(customMax[j]-maxima[i]);
      if(minDis > distance){
        minDis = distance;
        idxRef = j;
      }
    }
    maximaClustMap[i] = idxRef;
    //cout << "m i " << i << "  " << idxRef << "  " << maxima[i] << endl;

  }

  return labelCells(nearestMaximum, maximaClustMap, body);
}

// sets parameters that come from MGX
void RootCellProcessing::setParameter(bool outputType, int label1, int label2)
{

  (*config)[0].twoParts = outputType;

  //rootData.hasRadicle = outputType; 

  // get the two selected cells in order
  int idx1 = (*config)[0].labelIdxMap[label1];
  int idx2 = (*config)[0].labelIdxMap[label2];
  if((*config)[0].sigma1 == 0) (*config)[0].sigma1 = 4.0;
  if((*config)[0].sigma2 == 0) (*config)[0].sigma2 = 6.0;

  (*config)[0].labelLastRootCapCell = label1;
  (*config)[0].labelFirstCell = label2;

  if((*data)[idx1].coordLong < (*data)[idx2].coordLong){
    (*config)[0].labelLastRootCapCell = label2;
    (*config)[0].labelFirstCell = label1;
  }

  

  initHeatMapDataStructure(mainBody, (*config)[0].sigma1, (*config)[0].part1ChosenX, (*config)[0].part1ChosenY);
  if((*config)[0].twoParts)
    initHeatMapDataStructure(radicle, (*config)[0].sigma2, (*config)[0].part2ChosenX, (*config)[0].part2ChosenY);

}

// finds highest heatmap value for scaling the colours in the heatmap
// actually finds the 3 highest values and return the third value to avoid extreme outliers
double RootCellProcessing::findHighMax(heatMapDataStructure& body)
{

  std::map<int,Point2d> maxVec = body.maximaHeatMapAll;
  cellMorphLandscape landsc = body.heatMap;

  double high1 = 0, high2 = 0, high3 = 0;

  int maxVecSize = maxVec.size();

  for(int i=0; i<maxVecSize; i++)
    {
    Point2d currentPos = maxVec[i];
    int x = currentPos[0];
    int y = currentPos[1];
    double maxima = landsc[x][y];
      if(high1 < maxima)
      {
        high3 = high2; 
        high2 = high1;
        high1 = maxima;
      } else if(high2 < maxima)
      {
        high3 = high2; 
        high2 = maxima;
      } else if(high3 < maxima)
      {
        high3 = maxima; 
      }
    }
    //cout << "highs3  " << high1 << "  " << high2 << "  " << high3 << endl;
  return high3;

}

// sets the X coord of the heatmap
void RootCellProcessing::setHeatmap(heatMapDataStructure& body, QString optionGUI, int& chosenOption, bool x)
{

  QString option;

  if(optionGUI == "Radial Cell Length") {
    //body.x = rootData.lengthRad;
    chosenOption = 4;
    option = optionGUI;
  } else if(optionGUI == "Circumferential Cell Length") {
    //body.x = rootData.lengthCirc;
    chosenOption = 5;
    option = optionGUI;
  } else if(optionGUI == "Longitudinal Cell Length") {
    //body.x = rootData.lengthLong;
    chosenOption = 6;
    option = optionGUI;
  } else if(optionGUI == "Radial Coordinate") {
    //body.x = rootData.scaledRadialDis;
    chosenOption = 1;
    option = "Radial";
  } else if(optionGUI == "Circumferential Coordinate") {
    //body.x = rootData.azimuthal;
    chosenOption = 2;
    option = "Circumferential Angle";
  } else if(optionGUI == "Longitudinal Coordinate") {
    //body.x = rootData.arclengths;
    chosenOption = 3;
    option = "Longitudinal (Arclengths)";
  } else if(optionGUI == "Cell Volume") {
    //body.x = rootData.arclengths;
    chosenOption = 7;
    option = optionGUI;
  } else if(optionGUI == "Cell Wall Area") {
    //body.x = rootData.arclengths;
    chosenOption = 8;
    option = optionGUI;
  }

  if(x)
    body.x = getMapFromAttributes(option);
  else
    body.y = getMapFromAttributes(option);
}

// big procedure to generate the heatmap and label the cells
void RootCellProcessing::geometryBasedLabelling()
{
  // Function to label cell data based on steepest ascent on fuzzy histogram of
  // cell morphology data, matched to a learning set. The input plot_out gives
  // you the option to see the segmentation.

  // check where is the start and where the end of the root
  correctDirections();

  std::map<int,double> cellLabelMain, cellLabelRoot;

  //int numCells = (*data).size();

  if((*config)[0].twoParts){
  // embryo OR root cap
    radicle.sigma = (*config)[0].sigma2;
    radicle.heatMap = generateCellMorphologyLandscape(radicle);
    radicle.maximaHeatMapAll = findMaxima(radicle.heatMap);
    radicle.high = findHighMax(radicle);
    cellLabelRoot = clusterCells(radicle);

  }
  // main body for all root types
  mainBody.sigma = (*config)[0].sigma1;
  mainBody.heatMap = generateCellMorphologyLandscape(mainBody);
  mainBody.maximaHeatMapAll = findMaxima(mainBody.heatMap);
  mainBody.high = findHighMax(mainBody);
  cellLabelMain = clusterCells(mainBody);

  
  // now fill the data structure to save all labels
  for(int i = 0; i<numCells; i++){
    int currentLabel = labelMap[i];//(*data)[i].cellLabel;//rootData.uniqueLabels[i];
    if((mainBody.cells.find(currentLabel) != mainBody.cells.end()) and !mainBody.preselect[currentLabel]){
      (*data)[i].parentLabel = cellLabelMain[currentLabel]+1;
    } else if((radicle.cells.find(currentLabel) != radicle.cells.end()) and!radicle.preselect[currentLabel]){
      (*data)[i].parentLabel =  cellLabelRoot[currentLabel]+1;
    } else{
    }
  }

}

// reset preselection and clusters
void RootCellProcessing::resetHeat(heatMapDataStructure& body)
{
  std::cout << "Reset Preselection" << std::endl;
  for(int i = 0; i<rootData.numCells; i++){
    int currentLabel = rootData.uniqueLabels[i];
    if(body.cells.find(currentLabel) != body.cells.end()){
      body.preselect[currentLabel] = false;
    }
  }
  std::cout << "Reset Clusters" << std::endl;
  std::map<int,Point2d> newClusters;
  body.maximaHeatMapSelected = newClusters;
}

// take all currently selected clusters, look for nearest maximum and handle them as preselected and remove them from the heatmap
void RootCellProcessing::preselectCells(heatMapDataStructure& body)
{

  std::map<int,Point2d> maxPre = body.maximaHeatMapSelected;
  IntIntMap maxPreLabel = body.maximaHeatMapSelectedLabel;

  cellMorphLandscape landsc = body.heatMap;
  std::set<int> rootArea = body.cells;
  std::map<int,double> x = body.x;
  std::map<int,double> y = body.y;
  std::map<int,Point2d> maxima = body.maximaHeatMapAll;

  // relate each cell to the nearest maximum
  std::map<int,double> nearestMaximum = relateCellstoMaxima(body, landsc, x, y, maxima);

  // find the nearest maxima to the selected postions

  int numMaxPre = maxPre.size();
  int numMax = maxima.size();

  IntIntMap preToMax;

  for(int i = 0; i<numMaxPre; i++){
    double minDis = 1E20;
    int idxRef = 0;
    for(int j = 0; j<numMax; j++){
      double distance = norm(maxPre[i]-maxima[j]);
      if(minDis > distance){
        minDis = distance;
        idxRef = j;
      }
    }
    preToMax[i] = idxRef;
  }


  //int numCells = (*data).size();

  for(int i = 0; i<numCells; i++){
    int currentLabel = labelMap[i];//(*data)[i].cellLabel;//rootData.uniqueLabels[i];
    if(body.cells.find(currentLabel) != body.cells.end()){
      for(uint j = 0; j<preToMax.size(); j++){
        if(nearestMaximum[currentLabel] == preToMax[j]){
          (*data)[i].parentLabel = maxPreLabel[j]+1;
          body.preselect[currentLabel] = true;
        }
      }
    }
  }

  std::map<int,Point2d> newMax;
  body.maximaHeatMapSelected = newMax;


}

// label columella cells
// cells with small wall area to neighbors become root cap
// remaining cells: small become vasc, big col
void RootCellProcessing::assignColumella(int labelRootCap, int labelCol, double rat_val, double sca_val)//, int labelVasc)
{
  intToVec cellConnections;// = rcp.rootData.cellConnections;
  intToVec cellWallAreaPaired;// = rcp.rootData.cellWallAreas;

  //int numCells = (*data).size();

  // create a connection list
  forall(const IntIntDouPair &pr, (*config)[0].wallArea){
    cellConnections[pr.first.first].push_back(pr.first.second);
    cellWallAreaPaired[pr.first.first].push_back(pr.second);
  }

  double maxVol = 0;

  // calc max volume of the columella cells
  for(int i = 0; i<numCells; i++){
    //int currentLabel = labelMap[i];//(*data)[i].cellLabel;
    if((*data)[i].isColumella and (*data)[i].volume > maxVol)
      maxVol = (*data)[i].volume;
  }

  std::map<int,double> wallShared;
  std::map<int,double> ratio;

  // calculate ratio of total shared wall area to wall area
  for(int i = 0; i<numCells; i++){
    int currentLabel = labelMap[i];//(*data)[i].cellLabel;
    if(!(*data)[i].isColumella) continue;
    for(uint j = 0; j<cellConnections[currentLabel].size(); j++){
      wallShared[currentLabel]+=cellWallAreaPaired[currentLabel][j];
    }

    ratio[currentLabel] = wallShared[currentLabel] / (*data)[i].wallArea;
  }

  // find root cap cells
  std::map<int,bool> rootCapCells;

  // find vasculature cells
  //std::map<int,bool> vasc;
  //double volThreshold = 0.2;

  // all other cells
  std::map<int,bool> nonRootCap;

  for(int i = 0; i<numCells; i++){
    int currentLabel = labelMap[i];//(*data)[i].cellLabel;//rootData.uniqueLabels[i];
    if(!(*data)[i].isColumella) continue;
    if(ratio[currentLabel] < rat_val and (*data)[i].coordRad > sca_val){
    // if small wallarea to neighbors and radial distance big enough -> root cap
      rootCapCells[currentLabel] = true;
    } else {//if(!vasc[currentLabel] && !rootCapCells[currentLabel]){
      nonRootCap[currentLabel] = true;
    }
    // "last ditch" in the matlab files
    if(nonRootCap[currentLabel] and (*data)[i].coordRad > 0.6){
      rootCapCells[currentLabel] = true;
      nonRootCap[currentLabel] = false;
    }

    if(rootCapCells[currentLabel]){
      (*data)[i].parentLabel = labelRootCap;
    } else {
      (*data)[i].parentLabel = labelCol;
    }

    }

}

// assign cortical cells subroutine to create the correctedlengths and to dffind the min and max of the arclengths
void RootCellProcessing::assignCortCellsGetMinMax(int labelCort, double& minS, double& maxS)
{

  //int numCells = (*data).size();

  //int sizeBez = rootData.bez.size();
  int sizeBez = (*config)[0].bez.size();
  double lengthBez = 0;

  for(int i = 1; i<sizeBez; i++){
    lengthBez+=norm((*config)[0].bez[i]-(*config)[0].bez[i-1]);
  }
cout << "lbez " <<  lengthBez << endl;
  std::map<int,Point3d> longDir;

  std::map<int,double> correctedLengths;

  minS = 1E20;
  maxS = -1E20;

  for(int i = 0; i<numCells; i++){
    int currentLabel = labelMap[i];//(*data)[i].cellLabel; // rootData.uniqueLabels[i];
    if((*data)[i].parentLabel == labelCort){
      Point3d bezInterpPoint = (*config)[0].diffBezInterp[currentLabel];

      bezInterpPoint = -bezInterpPoint/norm(bezInterpPoint);
      correctedLengths[currentLabel] = -(*data)[i].sizeLong*((*data)[i].dirLong*bezInterpPoint);
      (*data)[i].correctedLengths = correctedLengths[currentLabel];
      //cout << "corr " <<  (*data)[i].correctedLengths << "/" << (*data)[i].sizeLong << "/" << (*data)[i].dirLong << "/" << bezInterpPoint << endl;
      if(minS > (*data)[i].coordLong) minS = (*data)[i].coordLong;
      if(maxS < (*data)[i].coordLong) maxS = (*data)[i].coordLong;

    }
    
  }

  //rootData.correctedLengths = correctedLengths;
}

// assign the cortical cell number to each cell in the root
void RootCellProcessing::assignCortCells(int labelCort, std::vector<double>& sSmooth, std::map<int, Point3d>& bMapS, int length)
{

  //int numCells = (*data).size();

  // now spline fit x = arclengths / y = correctedLengths of labelCort cells
  int cortCells = 0;
  for(int i = 0; i<numCells; i++){
    //int currentLabel = labelMap[i];//(*data)[i].cellLabel; // rootData.uniqueLabels[i];
    if((*data)[i].parentLabel == labelCort){
      cortCells++;
    }
  }  
  //cout << "spline" << endl;
  double x[cortCells], y[cortCells];
  cortCells = 0;
  for(int i = 0; i<numCells; i++){
    //int currentLabel = labelMap[i];//(*data)[i].cellLabel; // rootData.uniqueLabels[i];
    if((*data)[i].parentLabel == labelCort){
      x[cortCells] = (*data)[i].coordLong;//rootData.arclengths[currentLabel];
      y[cortCells] = (*data)[i].correctedLengths;// rootData.correctedLengths[currentLabel];
      cortCells++;
    }
  }  

  double eps = 0.000001;

  // TODO refactor, use a map or something else that doesnt require bubble sort

  // sort the x values from lowest to highest (bubble sort implementation)
  // x vector has to be increasing for the gsl spline
  for(int i = 0; i<cortCells; i++){
    for(int j = 0; j<cortCells; j++){
      if(i>j and x[i] < x[j]){
        double a = x[i];
        x[i] = x[j];
        x[j] = a;
        a = y[i];
        y[i] = y[j];
        y[j] = a;
      }
      // equal x values have to be eliminated // not sure if still necessary
      if(x[i] == x[j] and i!=j){
        if(i<j) x[i]-=eps;
        if(i>j) x[j]-=eps;
      }
    }
  }


  const size_t ncoeffs = 12; // number of fit coef
  const size_t nbreak = 10; // ncoeff-2
  gsl_bspline_workspace *bw;
  gsl_vector *B;
  gsl_vector *c, *w;
  gsl_vector *xGsl, *yGsl;
  gsl_matrix *X, *cov;
  gsl_multifit_linear_workspace *mw;
  double chisq;

  // allocate cubic bspline workspace (k = 4)
  bw = gsl_bspline_alloc(4, nbreak);
  B = gsl_vector_alloc(ncoeffs);

  xGsl = gsl_vector_alloc(cortCells);
  yGsl = gsl_vector_alloc(cortCells);
  X = gsl_matrix_alloc(cortCells, ncoeffs);
  c = gsl_vector_alloc(ncoeffs);
  w = gsl_vector_alloc(cortCells);
  cov = gsl_matrix_alloc(ncoeffs, ncoeffs);
  mw = gsl_multifit_linear_alloc(cortCells, ncoeffs);

  // this is the data to be fitted
  for (int i = 0; i < cortCells; i++)
    {
      gsl_vector_set(xGsl, i, x[i]);
      gsl_vector_set(yGsl, i, y[i]);
      gsl_vector_set(w, i, 1.0);
    }

  // uniform breakpoints
  gsl_bspline_knots_uniform(x[0], x[cortCells-1], bw);

  // construct the fit matrix X
  for (int i = 0; i < cortCells; i++)
    {
      double xi = gsl_vector_get(xGsl, i);

      // compute B_j(xi) for all j
      gsl_bspline_eval(xi, B, bw);

      // fill in row i of X
      for (uint j = 0; j < ncoeffs; j++)
        {
          double Bj = gsl_vector_get(B, j);
          gsl_matrix_set(X, i, j, Bj);
        }
    }

  // do the fit
  gsl_multifit_wlinear(X, w, yGsl, c, cov, &chisq, mw);

  double totalLength = 0;

  std::vector<double> sSmooth2;

  for(int i = 1; i<length; i++){
    totalLength+=norm(bMapS[i]-bMapS[i-1]);
  }

  double divider = sSmooth[length-1]-sSmooth[0];
  if(divider<0) divider*= -1.0;

  double sMin = 1E20;
  double sMax = -1E20;

  for(int i = 0; i<length; i++){
    sSmooth2.push_back(sSmooth[i]*totalLength/divider);
    if(sSmooth2[i]<sMin) sMin = sSmooth2[i];
    if(sSmooth2[i]>sMax) sMax = sSmooth2[i];
  } 

  double sTest = sMin;
  int rounds = 0;
  std::vector<double> sBin;

  // go through whole root and create bins of the size of the average cortical cell size
  while(sTest <= sMax){
    double splineVal, yerr;
    gsl_bspline_eval(sTest*divider/totalLength, B, bw);
    gsl_multifit_linear_est(B, c, cov, &splineVal, &yerr);
    sTest+=splineVal;
    sBin.push_back(sTest);
    rounds++;
  }

  // assign the right bin to each cell
  for(int i = 0; i<numCells; i++){
    //int currentLabel = labelMap[i];//(*data)[i].cellLabel;
    //double currentCompare = rootData.arclengths[currentLabel]*totalLength/divider;
    double currentCompare = (*data)[i].coordLong*totalLength/divider;
    for(int j = 2; j<rounds; j++){
      if(currentCompare<sBin[j] and currentCompare>=sBin[j-1]){
        (*data)[i].associatedCorticalCell = j-1;
      }
    }
    if(currentCompare>=sBin[rounds-1]){
      (*data)[i].associatedCorticalCell = rounds-1;
    }
    if(currentCompare<sBin[0]){
      (*data)[i].associatedCorticalCell = 0;
    }
  }

  // free gsl structures
  gsl_bspline_free(bw);
  gsl_vector_free(B);
  gsl_vector_free(xGsl);
  gsl_vector_free(yGsl);
  gsl_matrix_free(X);
  gsl_vector_free(c);
  gsl_vector_free(w);
  gsl_matrix_free(cov);
  gsl_multifit_linear_free(mw);

}


void RootCellProcessing::createHeatMapNoCellAtlas(AttrMap<int, double>& dataX, AttrMap<int, double>& dataY)
{

  mainBody.gridSize = 160;
  if(mainBody.sigma==0) mainBody.sigma = 5.0;
  numCells = 0;
  std::map<int,int> newLabelMap;

  typedef std::pair<int, double> IntDoubleP;
  forall(IntDoubleP p, dataX){
    mainBody.x[p.first] = p.second;

    newLabelMap[numCells] = p.first;
    numCells++;
    mainBody.cells.insert(p.first);
    //cout << "data read " << p.first << "/" << p.second << endl;
  }
  numCells = 0;
  forall(IntDoubleP p, dataY){
    mainBody.y[p.first] = p.second;
    numCells++;
  }
  labelMap = newLabelMap;

  //std::cout << "generateCellMorphologyLandscape" << std::endl;
  mainBody.heatMap = generateCellMorphologyLandscape(mainBody);
  //std::cout << "findMaxima" << std::endl;
  mainBody.maximaHeatMapAll = findMaxima(mainBody.heatMap);
 // std::cout << "findHighMax" << std::endl;
  mainBody.high = findHighMax(mainBody);

  std::map<int,double> cellLabelMain;
  cellLabelMain = clusterCells(mainBody);

  // now fill the data structure to save all labels
  for(int i = 0; i<numCells; i++){
    int currentLabel = labelMap[i];//(*data)[i].cellLabel;
    //int currentLabel = labelMap[i];//(*data)[i].cellLabel;//rootData.uniqueLabels[i];
    //if((mainBody.cells.find(currentLabel) != mainBody.cells.end()) and !mainBody.preselect[currentLabel]){
      parentsMap[currentLabel] = cellLabelMain[currentLabel]+1;

      //cout << "cell label/parent label  " << p.first << "/" << parentsMap[p.first] << endl;
    //}
  }

}

void RootCellProcessing::associateNearestCluster(heatMapDataStructure& body)
{
  IntIntMap maximaLabel = body.maximaHeatMapSelectedLabel;
  std::map<int,Point2d> maxima = body.maximaHeatMapSelected;

  std::map<int,double> cellLabel;
  // generate data for the graphs

  //int numCells = (*data).size();

  for(int i = 0; i<numCells; i++){
    int currentLabel = labelMap[i];//(*data)[i].cellLabel;
    if(!body.preselect[currentLabel] and body.cells.find(currentLabel) != body.cells.end()){

      double minDis = 1E20;
      int minLabel = -1;
      forall(IntInt p, maximaLabel){
        double dis = norm(body.pointsPos[currentLabel] - maxima[p.first]);
        if(dis < minDis){
          minDis = dis;
          minLabel = p.second;
        }
      }
      body.pointsLabel[body.pointsPos[currentLabel]] = minLabel+1;

    }
  }
}

/*
  Point2d minMaxFromMap(AttrMap<int,double>& inputMap, std::set<int>& keys)
  {
    Point2d res(0,0);

    double min = 1E20;
    double max = -1E20;


    forall(int k, keys){
      max = inputMap[k]>max? inputMap[k] : max;
      min = inputMap[k]<min? inputMap[k] : min;
    }

    double diff = max-min;

    if(min<1E20 and max>-1E20) res = Point2d(min-0.05*diff,max+0.05*diff);
    return res;

  }

  void ClusterMap::initDataArray()
  {
   std::vector<std::vector<double> > viewDataNew(gridSize);
   forall(std::vector<double>& v, viewDataNew){
     std::vector<double> vv(gridSize);
     v = vv;
   }
   viewData = viewDataNew;
  }



  void ClusterMap::createCellView()
  {
 cout << "cell view " << cellLabels.size() << "/" << cluster.size() << endl;
    gridSize = 480;
    gridFactor = 1;

    if(largeDots){
      gridSize = 160;
      gridFactor = 3;
    }
    initDataArray();
 

    // find measure1 and measure2 in idx map
    int idx1 = attrIdxMap[selectedX];
    int idx2 = attrIdxMap[selectedY];
    if(idx1==0 or idx2==0) return;

    // vector idx starts from 0
    idx1--;
    idx2--;

    // find min and max value in maps
    Point2d minMax1 = xMinMax;//minMaxFromMap(cellFeatures[idx1], cellLabels);
    Point2d minMax2 = yMinMax;//minMaxFromMap(cellFeatures[idx2], cellLabels);

    if(!customMinMax){
      minMax1 = minMaxFromMap(cellFeatures[idx1], cellLabels);
      minMax2 = minMaxFromMap(cellFeatures[idx2], cellLabels);
      xMinMax = minMax1;
      yMinMax = minMax2;
    }



    // bin the values

    // step size of the heatmap
    double step1 = (minMax1.y()-minMax1.x())/(double)gridSize;
    double step2 = (minMax2.y()-minMax2.x())/(double)gridSize;

    // reset viewdata array
    //for(int i = 0; i<gridSize; i++){
    //  for(int j = 0; j<gridSize; j++){
    //    viewData[i][j] = 0;
    //  }
    //}
 cout << "cell view " << endl;
    // fill array with new data
    forall(int label, cellLabels){
      int binX = interpolateArrayIndex(cellFeatures[idx1][label],minMax1.x(),minMax1.y(),gridSize-1);
      int binY = interpolateArrayIndex(cellFeatures[idx2][label],minMax2.x(),minMax2.y(),gridSize-1);

      viewData[binX][binY]++;
      //std::cout << "fill  " << binX <<"/" << binY << std::endl;
      imageLabelMap[std::make_pair(binX,binY)] = label;
      //std::cout << "f2  " << cellFeatures[idx1][label] <<"/" << minMax1.y() << "//" << minMax1.x() << std::endl;
    }
    if(!customMinMax) calcHeatMax();
  }

  void ClusterMap::createHeatMap2D() 
  {
    gridSize = 160;
    gridFactor = 3;
    initDataArray();

    // find measure1 and measure2 in idx map
    int idx1 = attrIdxMap[selectedX] -1;
    int idx2 = attrIdxMap[selectedY] -1;
    int idxHeat = attrIdxMap[selectedHeat] -1;
    if(idx1==-1 or idx2==-1) return;

    std::cout << "2d " << idx1 << "/" << idx2 << "//" << idxHeat << "///" << cellFeatures.size() << std::endl;


    // find min and max value in maps
    Point2d minMax1 = xMinMax;//minMaxFromMap(cellFeatures[idx1], cellLabels);
    Point2d minMax2 = yMinMax;//minMaxFromMap(cellFeatures[idx2], cellLabels);

    if(!customMinMax){
      minMax1 = minMaxFromMap(cellFeatures[idx1], cellLabels);
      minMax2 = minMaxFromMap(cellFeatures[idx2], cellLabels);
      xMinMax = minMax1;
      yMinMax = minMax2;
      //std::cout << "in" << std::endl;
    }


    // bin the values

    // step size of the heatmap
    double step1 = (minMax1.y()-minMax1.x())/(double)gridSize;
    double step2 = (minMax2.y()-minMax2.x())/(double)gridSize;

    //std::cout << "2d " << step1 << "/" << step2 << std::endl;

    // reset viewdata array
    //for(int i = 0; i<gridSize; i++){
    //  for(int j = 0; j<gridSize; j++){
    //    viewData[i][j] = 0;
    //  }
    //}

    //std::cout << "2d " << step1 << "/" << step2 << std::endl;

    // fill array with new data
    forall(int label, cellLabels){
      // calc center of gaussian
      int binX = interpolateArrayIndex(cellFeatures[idx1][label],minMax1.x(),minMax1.y(),gridSize-1);
      int binY = interpolateArrayIndex(cellFeatures[idx2][label],minMax2.x(),minMax2.y(),gridSize-1);

      // cut off radius around the center

// cut off radius around the center
      int radiusX = std::max(10.0,sigma*2);
      int radiusY = std::max(10.0,sigma*2);

      if(binX == 0 or binX == gridSize-1) radiusX = 0;
      if(binY == 0 or binY == gridSize-1) radiusY = 0;

      for(int i = -radiusX; i<=radiusX; i++){
        for(int j = -radiusY; j<=radiusY; j++){
          // calc gaussian
          if(binX+i>=0 and binX+i<gridSize and binY+j>=0 and binY+j<gridSize){
            double xValue = interpolateArrayIndex(cellFeatures[idx1][label]+i*step1,minMax1.x(),minMax1.y(),gridSize-1);
            double yValue = interpolateArrayIndex(cellFeatures[idx2][label]+j*step2,minMax2.x(),minMax2.y(),gridSize-1);
            if(idxHeat == -1)
              viewData[binX+i][binY+j] += gauss2D(xValue, yValue, binX, binY, sigma, sigma);
            else
              viewData[binX+i][binY+j] += cellFeatures[idxHeat][label] * gauss2D(xValue, yValue, binX, binY, sigma, sigma);
          }
        }
      }
    }
    if(!customMinMax) calcHeatMax();
    findMaximaHeatMap();
    relateCellsToMaxima();
    std::cout << "done createHeatMap2D" << std::endl;
  }

  void ClusterMap::calcHeatMax()
  {
    int size = gridSize;
    double maxValue = -HUGE_VAL;//body.high;
    double sum, count;
    std::map<int,Point2d> maxVec;// = body.maximaHeatMapAll;
  
    for(int i = 0; i < size; i++) {
      for(int j = 0; j < size; j++) {
        if(maxValue < viewData[i][j]) maxValue = viewData[i][j];
        if(viewData[i][j] > 0.1){
          sum+=viewData[i][j];
          count++;
        }
      }
    }
    heatMax = maxValue/5.;    
  }

void ClusterMap::setActiveMeasures(QString measure1, QString measure2, QString measureHeat)
  {

    std::cout << "set active " << std::endl;
    if(mode2D){
      Point2d xMinMaxOld = xMinMax;
      Point2d yMinMaxOld = yMinMax;
      forall(ManualCluster& c, cluster){
        if(selectedX != measure1){
          xMinMax = minMaxFromMap(cellFeatures[attrIdxMap[measure1]-1], cellLabels);
          double value1 = c.val[attrIdxMap[selectedX]-1];
          double posX  = realToImageCoord(xMinMaxOld, value1, gridSize);
          c.val[attrIdxMap[selectedX]-1] = 0;
          c.val[attrIdxMap[measure1]-1] = imageToRealCoord(xMinMax, posX, gridSize);
          //std::cout << "x new " << imageToRealCoord(xMinMax, posX, gridSize) << std::endl;
        }
        if(selectedY != measure2){
          yMinMax = minMaxFromMap(cellFeatures[attrIdxMap[measure2]-1], cellLabels);
          std::cout << "y " << std::endl;
          double value2 = c.val[attrIdxMap[selectedY]-1];
          double posY  = realToImageCoord(yMinMaxOld, value2, gridSize);
          c.val[attrIdxMap[selectedY]-1] = 0;
          c.val[attrIdxMap[measure2]-1] = imageToRealCoord(yMinMax, posY, gridSize);
        }
      }

    } else {
      std::cout << "nD " << std::endl;
    }

    selectedX = measure1;
    selectedY = measure2;
    selectedHeat = measureHeat;

  }

  QImage ClusterMap::createImage(bool heatMapMode){

    // init variables
    int size = gridSize;
    std::map<int,Point2d> maxVec;// = body.maximaHeatMapAll;

    // init image and colours
    int factor = gridFactor; // scaling factor
    QImage image( size*factor, size*factor, QImage::Format_ARGB32 );
    QRgb value;
    QRgb white = qRgba(255,255,255,255);
    QRgb grey1 = qRgba(200,200,200,255);
    QRgb purple = qRgba(120,0,120,255);
    QRgb black = qRgba(0,0,0,255);

    value = grey1;

    int yOffset = size*factor-1;

    // draw the heatmap and (if option is checked) small points for the cells
    for(int i = 0; i < size; i++) {
      for(int j = 0; j < size; j++) {
        //std::cout << "i/j/data " << i << "/" << j << "/" << data[i][j] << std::endl;
        for(int k = 0; k < factor; k++) {
          for(int l = 0; l < factor; l++) {
            if(heatMapMode) value = calcRGB(heatMax,viewData[i][j]);
            else {
              if(showParentLabels){
                value = getColorFromLabel(parentLabels[imageLabelMap[std::make_pair(i,j)]]);
              } else {
                value = calcRGB(heatMax,viewData[i][j]);
              }
              if(viewData[i][j] == 0) value = black;
            }
            image.setPixel(i*factor+k, yOffset-(j*factor+l), value);
          }
        }
      }
    }

    if(!heatMapMode) return image;

    // draw a small cross on all heatmap maxima
    int maxSize = maximaHeatMap.size();
    for(int i = 0; i < maxSize; i++) {
      Point2d currentPoint = maximaHeatMap[i];
      drawCross(image, currentPoint[0]*factor+factor/2, yOffset-currentPoint[1]*factor-factor/2, size*factor, 2, white);
    }

    return image;
  }

  double imageToRealCoord(Point2d minMax, double value, double gridSize)
  {
    double diff = minMax.y()-minMax.x();
    return value/gridSize * diff + minMax.x();

  }

  double realToImageCoord(Point2d minMax, double value, double gridSize)
  {
    double diff = minMax.y()-minMax.x();
    if(diff == 0) return 0;
    return ((value - minMax.x())/diff * gridSize);
  }

  void ClusterMap::addCluster(QString measure1, QString measure2, double value1, double value2)
  {
    ManualCluster newCluster;

    newCluster.val[attrIdxMap[measure1]-1] = imageToRealCoord(xMinMax, value1, gridSize);
    newCluster.val[attrIdxMap[measure2]-1] = imageToRealCoord(yMinMax, value2, gridSize);

    newCluster.label = 0;

    cluster.push_back(newCluster);
    std::cout << "new cluster " << value1 << "/" << value2 << std::endl;



    std::cout << "new cluster " << imageToRealCoord(xMinMax, value1, gridSize) << "/" << imageToRealCoord(yMinMax, value2, gridSize) << std::endl;


  }
  void ClusterMap::updateCluster(int clusterIdx, QString measure, double value, bool xCoord)
  {
    if(xCoord) cluster[clusterIdx].val[attrIdxMap[measure]-1] = imageToRealCoord(xMinMax, value, gridSize);
    else cluster[clusterIdx].val[attrIdxMap[measure]-1] = imageToRealCoord(yMinMax, value, gridSize);
  }

// checks whether there exists a maximum in maxVec close to the pos coordinates
// writes the result in maxIdx
bool ClusterMap::nearbyManualCluster(QString measure1, QString measure2, Point2d pos, int& maxIdx)
  {
    double minDis = 1E20;

    for(size_t i =0; i<cluster.size(); i++){
      ManualCluster& c = cluster[i];
    //forall(const ManualCluster& c, cluster){
      double realX = c.val[attrIdxMap[measure1]-1];
      double realY = c.val[attrIdxMap[measure2]-1];
      double imageX = interpolateArrayIndex(realX, xMinMax.x(), xMinMax.y(), gridSize-1);
      double imageY = interpolateArrayIndex(realY, yMinMax.x(), yMinMax.y(), gridSize-1);
      Point2d posImage(imageX,imageY);
      double dis = norm(pos - posImage);
      if(dis < minDis){
        minDis = dis;
        maxIdx = i;
      }
    }

    if(minDis < 10){
      return true;
    } else {
      maxIdx = -1;
      return false;
    }
  }

  // return the attr value for a specific label and measure
  double ClusterMap::getValue(QString measure, int label)
  {

    int idx = attrIdxMap[measure];
    if(idx==0) return 0.;
    idx--;
    return cellFeatures[idx][label];
  }

  double ClusterMap::getClusterValue(QString measure, int clusterIdx)
  {

    int idx = attrIdxMap[measure];
    if(idx==0) return 0.;
    idx--;
    return cluster[clusterIdx].val[idx];
  }

  void ClusterMap::setParentsCellView()
  {
    cout << "cell view " << cellLabels.size() << "/" << cluster.size() << endl;
    std::map<int,int> newParents;

    forall(int l, cellLabels){

    // find measure1 and measure2 in idx map
      int idx1 = attrIdxMap[selectedX];
      int idx2 = attrIdxMap[selectedY];
      if(idx1==0 or idx2==0) return;

      // vector idx starts from 0
      idx1--;
      idx2--;

//cout << "cell label " << l << "/" << attrIdxMap[measure1] << "/" << attrIdxMap[measure2] << "/" << cellFeatures.size() << endl;
      double realX = getValue(selectedX, l);//cellFeatures[idx1][l];
      double realY = getValue(selectedY, l);//cellFeatures[idx2][l];
//cout << "c " << realX << "/" << realY << endl;
      double imageX = realToImageCoord(xMinMax, realX, gridSize-1);
      double imageY = realToImageCoord(yMinMax, realY, gridSize-1);

      Point2d pCell(imageX,imageY);

      double minDis = HUGE_VAL;
      int minIdx = -1;
//cout << "c " << pCell << endl;
      for(size_t i =0; i<cluster.size(); i++){
        ManualCluster& c = cluster[i];
        double realCX = c.val[idx1];
        double realCY = c.val[idx2];
//cout << "cl " << realCX << "/" << realCY << endl;
        double imageCX = realToImageCoord(xMinMax, realCX, gridSize-1);
        double imageCY = realToImageCoord(yMinMax, realCY, gridSize-1);

        Point2d pCluster(imageCX,imageCY);
//cout << "cl " << pCluster << endl;
        double dis = norm(pCell - pCluster);
        if(dis<minDis){
          minDis = dis;
          minIdx = i;
        }
      }


      if(minIdx > -1)
        newParents[l] = cluster[minIdx].label+1;

        //cout << "parent " << minIdx << endl;

    }


    parentLabels = newParents;
  }

  void ClusterMap::setParentsHeatMapView()
  {
    
    std::map<int,int> newParents;

    // relate max in heatmap to manual clusters

    for(size_t i =0; i<maximaHeatMap.size(); i++){
      Point2d maxHM = maximaHeatMap[i];

      double minDis = HUGE_VAL;
      int minIdx = -1;


      for(size_t j =0; j<cluster.size(); j++){
        ManualCluster& c = cluster[j];
        Point2d clusterPos(getClusterValue(selectedX,j),getClusterValue(selectedY,j));

        double imageCX = realToImageCoord(xMinMax, clusterPos.x(), gridSize-1);
        double imageCY = realToImageCoord(yMinMax, clusterPos.y(), gridSize-1);
        Point2d clusterPos2(imageCX,imageCY);

        double dis = norm(clusterPos2 - maxHM);
        if(dis<minDis){
          minDis = dis;
          minIdx = j;
        }
      }
      maximaIdxManualClusterIdxMap[i]=minIdx;
      //std::cout << "cl " << i << "/" << minIdx <<  "/" << maximaHeatMap[minIdx] << std::endl;

    }


    // relate single cells to manual cluster
    forall(int l, cellLabels){
      //std::cout << "la " << l << "/" << cellLabelMaximaIdxMap[l] <<  "/" << maximaIdxManualClusterIdxMap[cellLabelMaximaIdxMap[l]] << "/" << cluster[maximaIdxManualClusterIdxMap[cellLabelMaximaIdxMap[l]]].label << std::endl;
      newParents[l] = cluster[maximaIdxManualClusterIdxMap[cellLabelMaximaIdxMap[l]]].label+1;
    }
    parentLabels = newParents;
  }


// find all local maxima in a cellMorphLandscape (=heatmap) and save their positions
void ClusterMap::findMaximaHeatMap()
{
  std::vector<Point2d> maxima;

  int size = gridSize;

  // check every point of the landscape
  for(int i=0; i<size; i++){
    for(int j=0; j<size; j++){
      bool isMax = true;
      for(int in=-1; in<=1; in++){
        for(int jn=-1; jn<=1; jn++){
          if(i+in >= 0 and i+in < size and j+jn >= 0 and j+jn < size and (jn!=0 or in!=0))
            if(viewData[i][j] <= viewData[i+in][j+jn]) isMax = false;
        }
      }
      if(isMax){
        Point2d maxPos(i,j);
        maxima.push_back(maxPos);
      }
    }
  }

  maximaHeatMap = maxima;
  }



// relate each cell in the root to a maxima in the heatmap by following the steepest gradient upwards until a maximum is reached
void ClusterMap::relateCellsToMaxima()
{
 
  std::map<int,double> nearestMaximum;
  //std::set<int> rootArea = body.cells;
  int numMax = maximaHeatMap.size();
  double minX = xMinMax.x();/// body.heatMapMinMax[0];
  double maxX = xMinMax.y();//body.heatMapMinMax[1];
  double minY = yMinMax.x();//body.heatMapMinMax[2];
  double maxY = yMinMax.y();//body.heatMapMinMax[3];

  //int numCells = (*data).size();

  //for(int i = 0; i<numCells; i++){
  forall(int currentLabel, cellLabels){
   // int currentLabel = labelMap[i];//(*data)[i].cellLabel;//rootData.uniqueLabels[i];
    //if(rootArea.find(currentLabel) != rootArea.end()){
      Point2d currentPos;
      double minDis = 1E20;
      int minNum = 0;

      double imageX = realToImageCoord(xMinMax, getValue(selectedX,currentLabel), gridSize-1);
      double imageY = realToImageCoord(yMinMax, getValue(selectedY,currentLabel), gridSize-1);

      currentPos[0] = imageX;//interpolateArrayIndex(,minX,maxX,gridSize-1);
      currentPos[1] = imageY;//interpolateArrayIndex(getValue(measure1,currentLabel),minY,maxY,gridSize-1);

      for(int j = 0; j<numMax; j++){ // find nearest maximum
        Point2d clusterPos = maximaHeatMap[j];

        if(minDis > norm(currentPos - clusterPos)){
          minDis = norm(currentPos - clusterPos);
          minNum = j;
        }
      }

      int stopCounter = 0;

      // if nearest maximum too far, then follow gradient
      while(minDis > sigma and stopCounter<100){
        stopCounter++;
        double maxValue = -HUGE_VAL;
        int changeX = 0;
        int changeY = 0;
        // find highest value in neighborhood and go there
        for(int i = -1; i<=1; i++){
          for(int j = -1; j<=1; j++){
            if(currentPos[0]+i >= 0 and currentPos[0]+i < gridSize and currentPos[1]+j >=0 and currentPos[1]+j < gridSize){
              if(maxValue < viewData[currentPos[0]+i][currentPos[1]+j]){
                maxValue = viewData[currentPos[0]+i][currentPos[1]+j];
                changeX = i;
                changeY = j;
              }
            }
          }
        }
        // find nearest max again
        minDis = HUGE_VAL;
        minNum = 0;
        currentPos[0] += changeX;
        currentPos[1] += changeY;
        for(int j = 0; j<numMax; j++){ // find nearest maximum
          Point2d clusterPos = maximaHeatMap[j];

          if(minDis > norm(currentPos - clusterPos)){
            minDis = norm(currentPos - clusterPos);
            minNum = j;
          }
        }
      }
      cellLabelMaximaIdxMap[currentLabel] = minNum;
    }
  //}

  //maximaIdxManualClusterIdxMap = cellLabelMaximaIdxMap;
  //return nearestMaximum;
}
*/

/*
void ClusterMap::convertClusters2DToND()
{
  mode2D = false;

  

}


// take current GUI position of manual clusters and save them to idx 0 and 1
void ClusterMap::convertClustersNDTo2D()
{
  resetCellFeatures();
  mode2D = true;

  // save current features
  AttrMap<int, double>& featureX = cellFeatures[attrIdxMap[selectedX]-1];
  AttrMap<int, double>& featureY = cellFeatures[attrIdxMap[selectedX]-1];

  addCellFeature(selectedX, featureX);
  addCellFeature(selectedY, featureY);
}
*/
}
