//
// This file is part of MorphoGraphX - http://www.MorphoGraphX.org
// Copyright (C) 2012-2016 Richard S. Smith and collaborators.
//
// If you use MorphoGraphX in your work, please cite:
//   http://dx.doi.org/10.7554/eLife.05864
//
// MorphoGraphX is free software, and is licensed under under the terms of the 
// GNU General (GPL) Public License version 2.0, http://www.gnu.org/licenses.
// 
#include <ImageData.hpp>

#include <ClipRegion.hpp>
#include <CutSurf.hpp>
#include <Dir.hpp>
#include <Geometry.hpp>
#include <Information.hpp>
#include <Insert.hpp>
#include <Mesh.hpp>
#include <Misc.hpp>
#include <Progress.hpp>
#include <Vector.hpp>
#include <ColorMap.hpp>
#include <TransferFunctionDlg.hpp>

#include <limits>
#include <math.h>

namespace mgx 
{
  typedef Vector<16, double> Point16d;
  
  using Information::err;
  
  using namespace std;
  
  // Static class varaibles
  int ImgData::VoxelEditRadius = 25;
  int ImgData::VoxelEditMaxPix = 1024000;
  int ImgData::ClearColor = 0;
  uint ImgData::Slices = 500u;
  uint ImgData::TileCount = 20u;
  Point3u ImgData::MaxTexSize = Point3u(2048u, 2048u, 2048u);
  float ImgData::SurfOffset = .01f;
  bool ImgData::SeedStack = false;
  
  float ImgData::DrawNormals = 0.0f;
  float ImgData::DrawZeroLabels = 0.0f;
  float ImgData::DrawNhbds = 0.0f;
  float ImgData::DrawOffset = 2.0f;
  bool ImgData::DeleteBadVertex = false;
  bool ImgData::FillWorkData = false;
  std::vector<Colorf> ImgData::LabelColors; // = std::vector<Point3f>(1);
  
  bool ImgData::MeshSelect;
  ScaleBar ImgData::scaleBar;
  Colorbar ImgData::colorBar;
  
  float ImgData::MeshPointSize = 1.0f;
  float ImgData::MeshLineWidth = 1.0f;
  
  #ifdef WIN32
  #  define FileDialogOptions QFileDialog::DontUseNativeDialog
  #else
  const int FileDialogOptions = 0;
  #endif
  
  ImgData::ImgData(int id, QWidget* _parent)
    : StackId(id), Main16Bit(false), Work16Bit(false), triangleColorEditDlg(0),
      newMainColorMap(true), newWorkColorMap(true), newSurfColorMap(true), newHeatColorMap(true),
      mainDataTexId(0), workDataTexId(0), dataTexColor(0), surfTexId(0), heatTexId(0), colorMapTexId(0),
      imgTexId(0), mcmapTexId(0), wcmapTexId(0), labelTexId(0), lineId(0), 
      selFboId(0), selFboColor(0), selFboDepth(0),
      posVAid(0), nrmlVAid(0), selVAid(0), texVAid(0), imgVAid(0), 
      triangleColorVAid(0), cellColorVAid(0),
      pntsVAid(0), pcolVAid(0), pselVAid(0), 
      cellGraphPosVAid(0), cellGraphColorVAid(0),
      axisPosVAid(0), axisColorVAid(0),
      texMap(0), parent(_parent), min(0), max(1),
      labelWallBordMin(0), labelWallBordMax(1), pixelRadius(1), marchLabel(0),
      meshShift(0), marchData(0), changed(false), pixelsChanged(false)
  {
    parent = _parent;
    surfTransferDlg = new TransferFunctionDlg(parent);
    heatTransferDlg = new TransferFunctionDlg(parent);
    workTransferDlg = new TransferFunctionDlg(parent);
    mainTransferDlg = new TransferFunctionDlg(parent);
  }
  
  void ImgData::init(Stack* s, Mesh* m)
  {
    stack = s;
    mesh = m;
    mainTransferDlg->setDefaultTransferFunction(stack->main()->transferFct());
    workTransferDlg->setDefaultTransferFunction(stack->work()->transferFct());
    surfTransferDlg->setDefaultTransferFunction(mesh->surfFct());
    heatTransferDlg->setDefaultTransferFunction(mesh->heatFct());
    updateMainColorMap();
    updateWorkColorMap();
    updateSurfColorMap();
    updateHeatColorMap();
  }
  
  ImgData::~ImgData()
  {
    // Delete VBOs
    if(posVAid)
      glfuncs->glDeleteBuffers(1, &posVAid);
    if(nrmlVAid)
      glfuncs->glDeleteBuffers(1, &nrmlVAid);
    if(selVAid)
      glfuncs->glDeleteBuffers(1, &selVAid);
    if(texVAid)
      glfuncs->glDeleteBuffers(1, &texVAid);
    if(imgVAid)
      glfuncs->glDeleteBuffers(1, &imgVAid);

    if(pntsVAid)
      glfuncs->glDeleteBuffers(1, &pntsVAid);
    if(pcolVAid)
      glfuncs->glDeleteBuffers(1, &pcolVAid);
    if(pselVAid)
      glfuncs->glDeleteBuffers(1, &pselVAid);

    if(cellGraphPosVAid)
      glfuncs->glDeleteBuffers(1, &cellGraphPosVAid);   
    if(cellGraphColorVAid)
      glfuncs->glDeleteBuffers(1, &cellGraphColorVAid); 

    if(axisPosVAid)
      glfuncs->glDeleteBuffers(1, &axisPosVAid); 
    if(axisColorVAid)
      glfuncs->glDeleteBuffers(1, &axisColorVAid); 
  
    // Delete textures
    if(mainDataTexId)
      glfuncs->glDeleteTextures(1, &mainDataTexId);
    if(workDataTexId)
      glfuncs->glDeleteTextures(1, &workDataTexId);
    if(surfTexId)
      glfuncs->glDeleteTextures(1, &surfTexId);
    if(heatTexId)
      glfuncs->glDeleteTextures(1, &heatTexId);
    if(colorMapTexId)
      glDeleteTextures(1, &colorMapTexId);
    if(imgTexId)
      glfuncs->glDeleteTextures(1, &imgTexId);
  
    delete mainTransferDlg;
    delete workTransferDlg;
    delete surfTransferDlg;
    delete heatTransferDlg;
  }
  
  // Set parameters for function
  void ImgData::readParms(Parms& parms, QString section)
  {
    Section = section;
  
    QString StackFile, WorkStackFile, MeshFile, ParentFile;
    parms(section, "StackFile", StackFile, QString(""));
    if(!StackFile.isEmpty())
      StackFile = absoluteFilePath(StackFile);
    stack->main()->setFile(StackFile);
    parms(section, "WorkStackFile", WorkStackFile, QString(""));
    if(!WorkStackFile.isEmpty())
      WorkStackFile = absoluteFilePath(WorkStackFile);
    stack->work()->setFile(WorkStackFile);
    parms(section, "MeshFile", MeshFile, QString(""));
    if(!MeshFile.isEmpty())
      MeshFile = absoluteFilePath(MeshFile);
    mesh->setFile(MeshFile);

    bool MainShow;
    parms(section, "MainShow", MainShow, true);
    if(MainShow)
      stack->main()->show();
    else
      stack->main()->hide();
    float MainBright;
    parms(section, "MainBright", MainBright, 1.0f);
    stack->main()->setBrightness(MainBright);
    float MainOpacity;
    parms(section, "MainOpacity", MainOpacity, 0.8f);
    stack->main()->setOpacity(MainOpacity);
    parms(section, "Main16Bit", Main16Bit, false);
    bool MainLabels;
    parms(section, "MainLabels", MainLabels, false);
    stack->main()->setLabels(MainLabels);
  
    bool WorkShow;
    parms(section, "WorkShow", WorkShow, false);
    if(WorkShow)
      stack->work()->show();
    else
      stack->work()->hide();
    float WorkBright;
    parms(section, "WorkBright", WorkBright, 1.0f);
    stack->work()->setBrightness(WorkBright);
    float WorkOpacity;
    parms(section, "WorkOpacity", WorkOpacity, 0.8f);
    stack->work()->setOpacity(WorkOpacity);
    parms(section, "Work16Bit", Work16Bit, false);
    bool WorkLabels;
    parms(section, "WorkLabels", WorkLabels, false);
    stack->work()->setLabels(WorkLabels);
  
    if(WorkStackFile == StackFile) {
      stack->work()->setFile(QString());
      if(WorkShow) {
        stack->main()->show();
        stack->work()->hide();
      }
    }
  
    bool SurfShow;
    parms(section, "SurfShow", SurfShow, true);
    mesh->setShowSurface(SurfShow);
    float SurfBright;
    parms(section, "SurfBright", SurfBright, 1.0f);
    mesh->setBrightness(SurfBright);
    float SurfOpacity;
    parms(section, "SurfOpacity", SurfOpacity, 1.0f);
    mesh->setOpacity(SurfOpacity);
    bool SurfCull, SurfBlend;
    parms(section, "SurfBlend", SurfBlend, false);
    parms(section, "SurfCull", SurfCull, true);
    mesh->setCulling(SurfCull);
    mesh->setBlending(SurfBlend);
  
    int SurfView;
    parms(section, "SurfView", SurfView, (int)Mesh::SURF_LABEL);
    if(SurfView >= Mesh::SURF_VIEW_COUNT or SurfView < 0)
      SurfView = (int)Mesh::SURF_LABEL;
    mesh->setSurfView((Mesh::SurfView)SurfView);

    QString view;
    parms(section, "SurfVertexView", view, QString("Signal"));
    mesh->setShowVertex(view, false);
    parms(section, "SurfTriangleView", view, QString("Triangle Value"));
    mesh->setShowTriangle(view, false);
    parms(section, "SurfLabelView", view, QString("Label"));
    mesh->setShowLabel(view, false);

		bool SurfParents;
    parms(section, "SurfParents", SurfParents, false);
    mesh->setUseParents(SurfParents);
  
    bool MeshShow, MeshPoints, MeshLines;
    parms(section, "MeshShow", MeshShow, true);
    mesh->setShowMesh(MeshShow);
    parms(section, "MeshPoints", MeshPoints, true);
    mesh->setShowMeshPoints(MeshPoints);
    parms(section, "MeshLines", MeshLines, true);
    mesh->setShowMeshLines(MeshLines);
    QString mvm;
    parms(section, "MeshViewMode", mvm, QString("Selected"));
    mesh->setMeshView(mvm, false);

    bool AxisShow;
    parms(section, "AxisShow", AxisShow, false);
    mesh->setShowAxis(AxisShow);
    QString avm;
    parms(section, "AxisViewMode", avm, QString("Cell Axis"));
    mesh->setAxisView(avm, false);
  
    bool CellMap, ShowTrans, ShowBBox;
    Point3f Scale = Point3f(1.0, 1.0, 1.0);
    // float Scale;
    parms(section, "CellMap", CellMap, false);
    mesh->setShowMeshCellMap(CellMap);
    parms(section, "ShowTrans", ShowTrans, false);
    stack->setShowTrans(ShowTrans);
    parms(section, "ShowBBox", ShowBBox, false);
    stack->setShowBBox(ShowBBox);
    // Show scaling is turned off by default
    parms(section, "ScaleX", Scale.x(), 1.0f);
    parms(section, "ScaleY", Scale.y(), 1.0f);
    parms(section, "ScaleZ", Scale.z(), 1.0f);

    bool LoadMeshScaleUm, LoadMeshTransform;
    parms(section, "LoadMeshScaleUm", LoadMeshScaleUm, true);
    mesh->setScaled(LoadMeshScaleUm);
    parms(section, "LoadMeshTransform", LoadMeshTransform, false);
    mesh->setTransformed(LoadMeshTransform);
    parms(section, "MeshShift", meshShift, 0.f);

    QString sf;
    parms(section, "MainColorMap", sf, QString());
    TransferFunction mainTransferFct = stack->main()->transferFct();
    mainTransferFct.clear();
    if(not sf.isEmpty()) {
      QString qsf = unshield(sf);
      mainTransferFct = TransferFunction::load(qsf);
    }
    if(mainTransferFct.empty()) {
      if(Section == "Stack1")
        mainTransferFct = TransferFunction::scale_green();
      else
        mainTransferFct = TransferFunction::scale_red();
    }
    mainTransferDlg->setTransferFunction(mainTransferFct);
    stack->main()->setTransferFct(mainTransferFct);
    updateMainColorMap();
  
    QString sf2;
    parms(section, "WorkColorMap", sf2, QString());
    TransferFunction workTransferFct = stack->work()->transferFct();
    workTransferFct.clear();
    if(not sf2.isEmpty()) {
      QString qsf = unshield(sf2);
      workTransferFct = TransferFunction::load(qsf);
    }
    if(workTransferFct.empty()) {
      if(Section == "Stack1")
        workTransferFct = TransferFunction::scale_cyan();
      else
        workTransferFct = TransferFunction::scale_yellow();
    }
    workTransferDlg->setTransferFunction(workTransferFct);
    stack->work()->setTransferFct(workTransferFct);
    updateWorkColorMap();
  
    QString sf3;
    parms(section, "SurfColorMap", sf3, QString());
    TransferFunction surfTransferFct = mesh->surfFct();
    surfTransferFct.clear();
    if(not sf3.isEmpty()) {
      QString qsf = unshield(sf3);
      surfTransferFct = TransferFunction::load(qsf);
    }
    if(surfTransferFct.empty())
      surfTransferFct = TransferFunction::scale_gray();
    surfTransferDlg->setTransferFunction(surfTransferFct);
    mesh->setSurfFct(surfTransferFct);
    updateSurfColorMap();
  
    QString sf4;
    parms(section, "HeatColorMap", sf4, QString());
    TransferFunction heatTransferFct = mesh->heatFct();
    heatTransferFct.clear();
    if(not sf4.isEmpty()) {
      QString qsf = unshield(sf4);
      heatTransferFct = TransferFunction::load(qsf);
    }
    if(heatTransferFct.empty())
      heatTransferFct = TransferFunction::jet();
    heatTransferDlg->setTransferFunction(heatTransferFct);
    mesh->setHeatFct(heatTransferFct);
    updateHeatColorMap();
  
    Matrix4d m = 1;
    parms(section, "Frame", m, m);
    stack->frame().setFromMatrix(m.c_data());
    m = 1;
    parms(section, "Trans", m, m);
    stack->trans().setFromMatrix(m.c_data());
  
    if(section == "Stack1") {
      MeshColor = Colors::Mesh1Color;
      MeshBorderColor = Colors::Mesh1BorderColor;
      // MeshPointColor = Colors::Mesh1PointColor;
      MeshSelectColor = Colors::Mesh1SelectColor;
    } else {
      MeshColor = Colors::Mesh2Color;
      MeshBorderColor = Colors::Mesh2BorderColor;
      // MeshPointColor = Colors::Mesh2PointColor;
      MeshSelectColor = Colors::Mesh2SelectColor;
    }
  }
  
  void ImgData::setSurfColorMap(const TransferFunction& fct) {
    surfTransferDlg->setTransferFunction(fct);
  }
  
  void ImgData::setHeatColorMap(const TransferFunction& fct) {
    heatTransferDlg->setTransferFunction(fct);
  }
  
  void ImgData::setWorkColorMap(const TransferFunction& fct) {
    workTransferDlg->setTransferFunction(fct);
  }
  
  void ImgData::setMainColorMap(const TransferFunction& fct) {
    mainTransferDlg->setTransferFunction(fct);
  }
  
  // Set parameters for function
  void ImgData::writeParms(QTextStream& pout, QString section)
  {
    Section = section;
  
    pout << endl;
    pout << "[" << section << "]" << endl;
  
    pout << "StackFile: " << stripCurrentDir(stack->main()->file()) << endl;
    pout << "WorkStackFile: " << stripCurrentDir(stack->work()->file()) << endl;
    pout << "MeshFile: " << stripCurrentDir(mesh->file()) << endl;
  
    pout << "MainShow: " << (stack->main()->isVisible() ? "true" : "false") << endl;
    pout << "MainBright:" << stack->main()->brightness() << endl;
    pout << "MainOpacity:" << stack->main()->opacity() << endl;
    pout << "Main16Bit:" << (Main16Bit ? "true" : "false") << endl;
    pout << "MainLabels:" << (stack->main()->labels() ? "true" : "false") << endl;
  
    pout << "WorkShow: " << (stack->work()->isVisible() ? "true" : "false") << endl;
    pout << "WorkBright:" << stack->work()->brightness() << endl;
    pout << "WorkOpacity:" << stack->work()->opacity() << endl;
    pout << "Work16Bit:" << (Work16Bit ? "true" : "false") << endl;
    pout << "WorkLabels:" << (stack->work()->labels() ? "true" : "false") << endl;
  
    pout << "SurfShow:" << (mesh->showSurface() ? "true" : "false") << endl;
    pout << "SurfBright:" << mesh->brightness() << endl;
    pout << "SurfOpacity:" << mesh->opacity() << endl;
    pout << "SurfBlend:" << (mesh->blending() ? "true" : "false") << endl;
    pout << "SurfCull:" << (mesh->culling() ? "true" : "false") << endl;
  
    pout << "SurfView:" << (int)mesh->surfView() << endl;
    pout << "SurfVertexView:" << mesh->surfVertexView() << endl;
    pout << "SurfTriangleView:" << mesh->surfTriangleView() << endl;
    pout << "SurfLabelView:" << mesh->surfLabelView() << endl;

    pout << "SurfParents:" << (mesh->useParents() ? "true" : "false") << endl;
  
    pout << "MeshShow:" << (mesh->showMesh() ? "true" : "false") << endl;
    pout << "MeshViewMode:" << mesh->meshView() << endl;
    pout << "AxisShow:" << (mesh->showAxis() ? "true" : "false") << endl;
    pout << "AxisViewMode:" << mesh->axisView() << endl;
    pout << "MeshPoints:" << (mesh->showMeshPoints() ? "true" : "false") << endl;
    pout << "MeshLines:" << (mesh->showMeshLines() ? "true" : "false") << endl;
  
    pout << "CellMap:" << (mesh->showMeshCellMap() ? "true" : "false") << endl;
    pout << "ShowTrans:" << (stack->showTrans() ? "true" : "false") << endl;
    pout << "ShowBBox:" << (stack->showBBox() ? "true" : "false") << endl;
    pout << "ShowScale:" << (stack->showScale() ? "true" : "false") << endl;
    pout << "TieScales:" << (stack->tieScales() ? "true" : "false") << endl;
    pout << "ScaleX:" << stack->scale().x() << endl;
    pout << "ScaleY:" << stack->scale().y() << endl;
    pout << "ScaleZ:" << stack->scale().z() << endl;
    pout << "Scale:" << norm(stack->scale()) << endl;

    pout << "LoadMeshScaleUm: " << (mesh->scaled() ? "true" : "false") << endl;
    pout << "LoadMeshTransform: " << (mesh->transformed() ? "true" : "false") << endl;
    pout << "MeshShift: " << meshShift << endl;  

    QString qsf = shield(stack->main()->transferFct().dump());
    pout << "MainColorMap: " << qsf << endl;
    qsf = shield(stack->work()->transferFct().dump());
    pout << "WorkColorMap: " << qsf << endl;
    qsf = shield(mesh->surfFct().dump());
    pout << "SurfColorMap: " << qsf << endl;
    qsf = shield(mesh->heatFct().dump());
    pout << "HeatColorMap: " << qsf << endl;
  
    pout << "Frame: " << Point16d(stack->frame().matrix()) << endl;
    pout << "Trans: " << Point16d(stack->trans().matrix()) << endl;
    pout << endl;
  }
  
  // Set stack sizes
  void ImgData::updateStackSize()
  {
    if(stack->empty())
      unloadTex();
  
    emit changeSize(stack->size(), stack->step(), stack->origin());
  }
  
  static void drawPlane(Shader* shader, Point3f imageSize, Point3u texSize, Point3f Origin, float scaling, int Slices,
                        GLboolean hasClip1, GLboolean hasClip2, GLboolean hasClip3, Point3f c1, Point3f c2, Point3f c3,
                        Point3f c4)
  {
    shader->setUniform("origin", GLSLValue(vec3(Origin)));
    shader->setUniform("imageSize", GLSLValue(vec3(imageSize)));
    shader->setUniform("scaling", GLSLValue(scaling));
    shader->setUniform("texSize", GLSLValue(ivec3(texSize)));
    shader->setUniform("slices", GLSLValue(float(Slices) / 10000.f));
    shader->setUniform("has_clipped", GLSLValue(ivec3(hasClip1, hasClip2, hasClip3)));
    shader->useShaders();
    shader->setupUniforms();
  
    glBegin(GL_QUADS);
    glColor4f(1, 1, 1, 1);
    glVertex3fv(c4.c_data());
    glVertex3fv(c3.c_data());
    glVertex3fv(c2.c_data());
    glVertex3fv(c1.c_data());
    glEnd();
  
    shader->stopUsingShaders();
  }
  
  void ImgData::setupVolumeShader(Shader& shader, int pos)
  {
    QString color_fct = "vec4 color(vec3 texCoord) { \n";
    if(isMainVisible() and isWorkVisible()) {
      if(stack->main()->labels())
        color_fct += "  vec4 c1 = (brightness > 0.0) ? premulColor(index_color1(texCoord)) : vec4(0,0,0,1);\n";
      else
        color_fct += "  vec4 c1 = (brightness > 0.0) ? premulColor(colormap_color1(texCoord)) : vec4(0,0,0,1);\n";
      if(stack->work()->labels())
        color_fct
          += "  vec4 c2 = (second_brightness > 0.0) ? premulColor(index_color2(texCoord)) : vec4(0,0,0,1);\n";
      else
        color_fct
          += "  vec4 c2 = (second_brightness > 0.0) ? premulColor(colormap_color2(texCoord)) : vec4(0,0,0,1);\n";
      color_fct += "  return mixColors(c1, c2);\n}\n";
    } else if(isMainVisible()) {
      if(stack->main()->labels())
        color_fct += "  return (brightness > 0.0) ? premulColor(index_color1(texCoord)) : vec4(0,0,0,1);\n}\n";
      else
        color_fct += "  return (brightness > 0.0) ? premulColor(colormap_color1(texCoord)) : vec4(0,0,0,1);\n}\n";
    } else if(isWorkVisible()) {
      if(stack->work()->labels())
        color_fct
          += "  return (second_brightness > 0.0) ? premulColor(index_color2(texCoord)) : vec4(0,0,0,1);\n}\n";
      else
        color_fct
          += "  return (second_brightness > 0.0) ? premulColor(colormap_color2(texCoord)) : vec4(0,0,0,1);\n}\n";
    } else
      color_fct += "  return vec4(0,0,0,1);\n}\n";
    shader.changeFragmentShaderCode(pos, color_fct);
  }
  
  void ImgData::setup3DRenderingData(Shader* shader)
  {
    // When drawing 3D seeds, we need to also draw the main stack
    if(stack->main()->isVisible() and mainDataTexId) {
      bind3DTex(mainDataTexId);
  
      setupMainColorMap();
  
      float brightness = powf(stack->main()->brightness(), 2.0f);
      float opacity = powf(stack->main()->opacity(), 2.0f);
      shader->setUniform("brightness", GLSLValue(brightness));
      shader->setUniform("opacity", GLSLValue(opacity));
    }
  
    if(stack->work()->isVisible() and workDataTexId) {
      bind3DTex(workDataTexId, Shader::AT_SECOND_TEX3D);
  
      setupWorkColorMap2();
  
      float brightness = powf(stack->work()->brightness(), 2.0f);
      float opacity = powf(stack->work()->opacity(), 2.0f);
      shader->setUniform("second_brightness", GLSLValue(brightness));
      shader->setUniform("second_opacity", GLSLValue(opacity));
    }
  }
  
  void ImgData::drawStack(Shader* shader)
  {
    if(!isVisible() or !(mainDataTexId or workDataTexId)
       or ((stack->main()->opacity() == 0 or !isMainVisible())and (stack->work()->opacity() == 0 or !isWorkVisible())))
      return;
  
    glEnable(GL_AUTO_NORMAL);
    glDisable(GL_LIGHTING);
    glDisable(GL_BLEND);
  
    glFrontFace(GL_CCW);
    glPolygonMode(GL_FRONT, GL_FILL);
  
    glDisable(GL_TEXTURE_1D);
    glDisable(GL_TEXTURE_2D);
    glDisable(GL_TEXTURE_3D);
  
    // Choose blending function
    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
    // glBlendFunc(GL_SRC_ALPHA, GL_ONE);
  
    // glBlendEquation(GL_MAX);
    // glBlendEquationSeparate(GL_MAX, GL_MIN);
  
    // Set matrix
    glMatrixMode(GL_MODELVIEW);
    glPushMatrix();
    glMultMatrixd(getFrame().worldMatrix());
  
    GLdouble mv[16], proj[16];
    GLint viewport[4];
    glGetDoublev(GL_MODELVIEW_MATRIX, mv);
    glGetDoublev(GL_PROJECTION_MATRIX, proj);
    glGetIntegerv(GL_VIEWPORT, viewport);
  
    // Get unit vectors in x,y,z direction
    Point3f x(mv[0], mv[4], mv[8]);
    Point3f y(mv[1], mv[5], mv[9]);
    Point3f z(mv[2], mv[6], mv[10]);
    // Vertex array for quads
    std::vector<Point3f> vtxVA;
  
    // Render object aligned or view aligned
    const Point3f& Size = stack->worldSize();
    double x2 = Size.x() * Size.x() / 4;
    double y2 = Size.y() * Size.y() / 4;
    double z2 = Size.z() * Size.z() / 4;
    float dis = sqrt(x2 + y2 + z2);
    z *= dis;
    x *= dis;
    y *= dis;
  
    Shader::activeTexture(Shader::AT_LABEL_TEX);
    glBindTexture(GL_TEXTURE_1D, labelTexId);
    Shader::activeTexture(Shader::AT_NONE);
  
    GLdouble winX = viewport[2];
    GLdouble winY = viewport[3];
  
    typedef Vector<3, GLdouble> Point3GLd;
    Point3GLd c1, c2, c3, c4;
    gluUnProject(0, 0, 0, mv, proj, viewport, &c1.x(), &c1.y(), &c1.z());
    gluUnProject(0, winY, 0, mv, proj, viewport, &c2.x(), &c2.y(), &c2.z());
    gluUnProject(winX, winY, 0, mv, proj, viewport, &c3.x(), &c3.y(), &c3.z());
    gluUnProject(winX, 0, 0, mv, proj, viewport, &c4.x(), &c4.y(), &c4.z());
  
    GLboolean clip1, clip2, clip3;
    glGetBooleanv(GL_CLIP_PLANE0, &clip1);
    glGetBooleanv(GL_CLIP_PLANE2, &clip2);
    glGetBooleanv(GL_CLIP_PLANE4, &clip3);
  
    setup3DRenderingData(shader);
  
    drawPlane(shader, Size, stack->size(), stack->origin(), stack->scale().x(),   // ALR: not sure what to do here
              Slices, clip1, clip2, clip3, Point3f(c1), Point3f(c2), Point3f(c3), Point3f(c4));
  
    // drawPlane(shader, Size, stack->size(), stack->origin(), stack->scale(),
    //          Slices,
    //          clip1, clip2, clip3,
    //          Point3f(c1), Point3f(c2), Point3f(c3), Point3f(c4));
    unbind3DTex();
  
    // glBlendEquation(GL_FUNC_ADD);
    glDisable(GL_BLEND);
    glMatrixMode(GL_MODELVIEW);
    glPopMatrix();
    REPORT_GL_ERROR("drawStack");
  }
  
  void ImgData::drawBBox()
  {
    if(!stack->showBBox() or !(mainDataTexId or workDataTexId))
      return;
    glDisable(GL_LIGHTING);
    glDisable(GL_TEXTURE_1D);
    glDisable(GL_TEXTURE_2D);
    glDisable(GL_TEXTURE_3D);
    glDisable(GL_BLEND);
    glDisable(GL_ALPHA_TEST);
    glDisable(GL_AUTO_NORMAL);
    glEnable(GL_DEPTH_TEST);
  
    glLineWidth(MeshLineWidth);
  
    glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
  
    glMatrixMode(GL_MODELVIEW);
    glPushMatrix();
    glMultMatrixd(getFrame().worldMatrix());
  
    Colorf c = Colors::getColor((StackId == 0) ? Colors::Stack1BBoxColor : Colors::Stack2BBoxColor);
  
    glColor4fv(c.c_data());
  
    const Point3f& Origin = stack->origin();
    Point3f c1 = Origin;
    Point3f c2 = Origin + stack->worldSize();
  
    glBegin(GL_QUADS);
    glVertex3f(c1.x(), c1.y(), c1.z());
    glVertex3f(c1.x(), c2.y(), c1.z());
    glVertex3f(c2.x(), c2.y(), c1.z());
    glVertex3f(c2.x(), c1.y(), c1.z());
  
    glVertex3f(c1.x(), c1.y(), c2.z());
    glVertex3f(c1.x(), c2.y(), c2.z());
    glVertex3f(c2.x(), c2.y(), c2.z());
    glVertex3f(c2.x(), c1.y(), c2.z());
    glEnd();
  
    glBegin(GL_LINES);
    glVertex3f(c1.x(), c1.y(), c1.z());
    glVertex3f(c1.x(), c1.y(), c2.z());
  
    glVertex3f(c1.x(), c2.y(), c1.z());
    glVertex3f(c1.x(), c2.y(), c2.z());
  
    glVertex3f(c2.x(), c2.y(), c1.z());
    glVertex3f(c2.x(), c2.y(), c2.z());
  
    glVertex3f(c2.x(), c1.y(), c1.z());
    glVertex3f(c2.x(), c1.y(), c2.z());
    glEnd();
  
    glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
    glMatrixMode(GL_MODELVIEW);
    glPopMatrix();
    REPORT_GL_ERROR("drawBBox");
  }
  
  // Draw mesh wireframe
  void ImgData::drawMesh()
  {
    if(!mesh->showMesh())
      return;
    if(pntsVA.size() == 0)
      return;
  
    glDisable(GL_LIGHTING);
    glDisable(GL_TEXTURE_1D);
    glDisable(GL_TEXTURE_2D);
    glDisable(GL_TEXTURE_3D);
    glDisable(GL_BLEND);
    glDisable(GL_ALPHA_TEST);
    glDisable(GL_AUTO_NORMAL);
    glEnable(GL_DEPTH_TEST);
    // glEnable(GL_MULTISAMPLE);
    glFrontFace(GL_CCW);
  
    glMatrixMode(GL_MODELVIEW);
    glPushMatrix();
    glMultMatrixd(getFrame().worldMatrix());
    // Scale stack if required
    if(stack->scale() != Point3f(1.0, 1.0, 1.0))
      glScaled(stack->scale().x(), stack->scale().y(), stack->scale().z());
    if(mesh->showMeshLines()) {
      glLineWidth(MeshLineWidth);
      glEnableClientState(GL_VERTEX_ARRAY);
      glEnableClientState(GL_COLOR_ARRAY); 
      if(mesh->meshView() == "Cell Graph"){
        glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, cellGraphPosVAid);
        glVertexPointer(3, GL_FLOAT, 0, 0);
        glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, cellGraphColorVAid);
        glColorPointer(3, GL_UNSIGNED_BYTE, 0, 0);
        glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);
        glDrawElements(GL_LINES, cellGraphLineVA.size(), GL_UNSIGNED_INT, &cellGraphLineVA[0]);
      } else {
        glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, pntsVAid);
        glVertexPointer(3, GL_FLOAT, 0, 0);
        glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, pcolVAid);
        glColorPointer(3, GL_UNSIGNED_BYTE, 0, 0);
        glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);
        if(mesh->meshView() == "All")
          glDrawElements(GL_LINES, lineVA.size(), GL_UNSIGNED_INT, &lineVA[0]);
        else if(mesh->meshView() == "Border")
          glDrawElements(GL_LINES, lbrdVA.size(), GL_UNSIGNED_INT, &lbrdVA[0]);
        else if(mesh->meshView() == "Cells")
          glDrawElements(GL_LINES, lcellVA.size(), GL_UNSIGNED_INT, &lcellVA[0]);
        else if(mesh->meshView() == "Selected")
          glDrawElements(GL_LINES, lselVA.size(), GL_UNSIGNED_INT, &lselVA[0]);
      }
      glDisableClientState(GL_VERTEX_ARRAY);
      glDisableClientState(GL_COLOR_ARRAY);
    }
  
    if(mesh->showMeshPoints()) {
      glEnableClientState(GL_VERTEX_ARRAY);
      glEnableClientState(GL_COLOR_ARRAY);
  
      glPointSize(MeshPointSize);
      if(mesh->meshView() == "Cell Graph"){
        glEnableClientState(GL_VERTEX_ARRAY);
        glEnableClientState(GL_COLOR_ARRAY);  
        glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, cellGraphPosVAid);
        glVertexPointer(3, GL_FLOAT, 0, 0);
        glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, cellGraphColorVAid);
        glColorPointer(3, GL_UNSIGNED_BYTE, 0, 0);
        glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);
        glDrawElements(GL_POINTS, cellGraphPointVA.size(), GL_UNSIGNED_INT, &cellGraphPointVA[0]);
      }
      if(mesh->meshView() == "All") {
        glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, pntsVAid);
        glVertexPointer(3, GL_FLOAT, 0, 0);
        glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, pcolVAid);
        glColorPointer(3, GL_UNSIGNED_BYTE, 0, 0);
        glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);
        glDrawArrays(GL_POINTS, 0, pntsVA.size());
      } else {
        glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, pntsVAid);
        glVertexPointer(3, GL_FLOAT, 0, 0);
        glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, pcolVAid);
        glColorPointer(3, GL_UNSIGNED_BYTE, 0, 0);
        glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);
        if(mesh->meshView() == "Border")
          glDrawElements(GL_POINTS, pbrdVA.size(), GL_UNSIGNED_INT, &pbrdVA[0]);
        else if(mesh->meshView() == "Cells")
          glDrawElements(GL_POINTS, pcellVA.size(), GL_UNSIGNED_INT, &pcellVA[0]);
        else if(mesh->meshView() == "Selected")
          glDrawElements(GL_POINTS, pselVA.size(), GL_UNSIGNED_INT, &pselVA[0]);
      }
      glDisableClientState(GL_VERTEX_ARRAY);
      glDisableClientState(GL_COLOR_ARRAY);
    }

    if(DrawNormals > 0) {
      glBegin(GL_LINES);
      for(uint i = 0; i < posVA.size(); i++) {
        glVertex3fv(posVA[i].data());
        glVertex3fv((posVA[i] + nrmlVA[i] * DrawNormals).data());
      }
      glEnd();
    }

// Not safe since mesh may be changing
//
//    if(DrawZeroLabels > 0) {
//      glLineWidth(MeshLineWidth);
//      Color3f col = Colors::getColor(MeshBorderColor);
//      glColor3fv(col.data());
//      vvGraph &S = mesh->graph();
//      glBegin(GL_LINES);
//      forall(const vertex& v, S) {
//        if(v->label != 0)
//          continue;
//        glVertex3fv(Point3f(v->pos).c_data());
//        glVertex3fv(Point3f(v->pos + DrawZeroLabels * v->nrml).c_data());
//      }
//      glEnd();
//    }
//
//    if(DrawNhbds > 0) {
//      glLineWidth(MeshLineWidth);
//      Color3f col = Colors::getColor(MeshBorderColor);
//      glColor3fv(col.data());
//      vvGraph &S = mesh->graph();
//      glBegin(GL_LINES);
//      forall(const vertex& v, S) {
//        float d = .1f;
//        forall(const vertex& n, S.neighbors(v)) {
//          Point3f a = (n->pos - v->pos).normalize();
//          glVertex3fv(Point3f(v->pos).c_data());
//          glVertex3fv(Point3f(v->pos + DrawNhbds * (d * a + .1f * v->nrml)).c_data());
//          d += .1f;
//        }
//      }
//      glEnd();
//    }

    /*
     * if(MeshBorderV) {
     *  glLineWidth(MeshLineWidth);
     *  Color3f col = Colors::getColor(MeshBorderColor);
     *  col /= 255.0;
     *  glColor3fv(col.data());
     *  glBegin(GL_LINES);
     *  forall(const vertex &v, S) {
     *    if(v->minb != 0) {
     *      vertex minb(v->minb);
     *      glVertex3fv(v->pos.c_data());
     *      glVertex3fv(minb->pos.c_data());
     *    }
     *  }
     *  glEnd();
     *}
     */
    glMatrixMode(GL_MODELVIEW);
    glPopMatrix();
    glDisable(GL_POLYGON_OFFSET_FILL);
    glDisable(GL_POLYGON_OFFSET_LINE);
    REPORT_GL_ERROR("drawMesh");
  }
  
  // Draw line in between pairs of junctions of two corresponding cell meshes (used in MeshProcessGrowth)
  void ImgData::drawVertexVertexLine(ImgData* otherstack)
  {
    const vvGraph& S1 = this->mesh->graph();
    const vvGraph& S2 = otherstack->mesh->graph();
  
    IntMatrix3fMap& VVCorrespondence = this->mesh->vvCorrespondence();
  
    if(!this->mesh->showMesh() and !otherstack->mesh->showMesh() and !this->mesh->showSurface()
       and !otherstack->mesh->showSurface())
      return;
  
    if(S1.empty() or S2.empty() or VVCorrespondence.empty()) {
      return;
    }
  
    bool showlines = this->mesh->showVVCorrespondence();
    if(!showlines)
      return;
  
    // stack matrices
    Matrix4d st1 = ~Matrix4d(getFrame().worldMatrix());
    Matrix4d st2 = ~Matrix4d(otherstack->getFrame().worldMatrix());
  
    glLineWidth(1);
    glColor3f(1., 1., 1.);
    glDisable(GL_LIGHTING);
    glBegin(GL_LINES);
    for(int i = 1; i <= (int)VVCorrespondence.size(); i++) {
      if(VVCorrespondence.count(i) > 0) {
        Point4d v1pos(VVCorrespondence[i][0]);
        Point4d v2pos(VVCorrespondence[i][1]);
        v1pos.t() = v2pos.t() = 1.0f;
        // Clipping planes, point normal form
        Point4f v1pos_transf(st1 * v1pos);
        Point4f v2pos_transf(st2 * v2pos);
        glVertex4fv(v1pos_transf.c_data());
        glVertex4fv(v2pos_transf.c_data());
      }
    }
    glEnd();
    glEnable(GL_LIGHTING);
    REPORT_GL_ERROR("drawVertexVertexLine");
  }
  
  // Draw cell axis (3 lines). The color is defined separately for each axis in each cell
  void ImgData::drawAxis()
  {
//    updateAxis(); // added to make the dropdown menu changer operational (StackxAxisView) // ALR: causes mesh drawing to be slow!

    if(!mesh->showAxis() or axisPosVA.size() == 0)
      return;
  
    const vvGraph& S = mesh->graph();
    if(S.empty()) 
      return;

    // Draw the lines
    REPORT_GL_ERROR("drawAxis");
  
    // to follow camera and object transformation:
    glMatrixMode(GL_MODELVIEW);
    glPushMatrix();
    glMultMatrixd(getFrame().worldMatrix());
    if(norm(stack->scale()) != 1.0)
      glScaled(stack->scale().x(), stack->scale().y(), stack->scale().z());
  
		glLineWidth(trim(mesh->axisWidth(), 1.0f, std::numeric_limits<float>::max()));
    REPORT_GL_ERROR("drawAxis");
    glDisable(GL_LIGHTING);
  
    glEnableClientState(GL_VERTEX_ARRAY);
    glEnableClientState(GL_COLOR_ARRAY);

    REPORT_GL_ERROR("drawAxis");
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, axisPosVAid);
    glVertexPointer(3, GL_FLOAT, 0, 0);
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, axisColorVAid);
    glColorPointer(3, GL_UNSIGNED_BYTE, 0, 0);
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);
    glDrawArrays(GL_LINES, 0, axisPosVA.size()); 

    REPORT_GL_ERROR("drawAxis");
    glDisableClientState(GL_VERTEX_ARRAY);
    glDisableClientState(GL_COLOR_ARRAY);
  
    glEnable(GL_LIGHTING);
  
    glPopMatrix();   // to follow camera and object transformation

    REPORT_GL_ERROR("drawAxis");
  }
  
  // Draw mesh shaded, with and without select
  void ImgData::drawSurf(bool select, Shader* textureShader, Shader* labelShader, 
                                          Shader *colorShader, Shader* volumeShader)
  {
    if((!select and !mesh->showSurface())
       or (select and MeshSelect))   // or (opaque_only and SurfBlend))
      return;

    if(idVA.size() == 0)
      return;

    // Are we drawing triangle colors
    bool drawTriColor = mesh->surfView() == Mesh::SURF_TRIANGLE and triangleColorVA.size() > 0;
    if(drawTriColor and triangleColorVA.size() != idVA.size())
      return;

    // Are we drawing cell colors
    bool drawCellColor = mesh->surfLabelView() == "Cell Color" 
                                     and mesh->surfView() == Mesh::SURF_LABEL; 
    if(drawCellColor and cellColorVA.size() != idVA.size())
      return;

    // Draw 2D image texture (like for Keyence)
    bool drawImageTex = mesh->surfVertexView() == "Image Texture";
    if(drawImageTex and !mesh->hasImgTex())
      return;

    // Color with the texture from the stack, in this case positions are texture coords
    bool drawStackTex = mesh->surfVertexView() == "Stack Texture";
    if(drawStackTex and stack->empty())
      return;

    // Bind select frame buffer
    if(select)
      startSelectFbo();

    glDisable(GL_LIGHTING);
    glDisable(GL_TEXTURE_1D);
    glDisable(GL_TEXTURE_2D);
    glDisable(GL_TEXTURE_3D);
    glDisable(GL_BLEND);
    glDisable(GL_AUTO_NORMAL);
    glDisable(GL_CULL_FACE);
    glEnable(GL_DEPTH_TEST);
    glFrontFace(GL_CCW);
    glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
  
    glEnableClientState(GL_VERTEX_ARRAY);

    if(select)
      glEnableClientState(GL_COLOR_ARRAY);
    else {
      glEnableClientState(GL_NORMAL_ARRAY);
      if(drawTriColor or drawCellColor)
        glEnableClientState(GL_COLOR_ARRAY);
      else
        glEnableClientState(GL_TEXTURE_COORD_ARRAY);
    }
  
    // First bind positions
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, posVAid);
    glVertexPointer(3, GL_FLOAT, 0, 0);
  
    Shader* shader = 0;
    Shader::activeTexture(Shader::AT_NONE);
  
    // Set matrix, scale if required
    glMatrixMode(GL_MODELVIEW);
    glPushMatrix();
    glMultMatrixd(getFrame().worldMatrix());
    // Scale stack if required
    if(stack->scale() != Point3f(1.0, 1.0, 1.0))
      glScaled(stack->scale().x(), stack->scale().y(), stack->scale().z());
  
    // Now do colors and/or texture
    if(select) {
      glClearColor(0, 0, 0, 1);
      glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
      glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, selVAid);
      glColorPointer(3, GL_UNSIGNED_BYTE, 0, 0);
    } else {
      glEnable(GL_LIGHTING);

      if(mesh->culling())
        glEnable(GL_CULL_FACE);

      glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, nrmlVAid);
      glNormalPointer(GL_FLOAT, 0, 0);

      if(drawTriColor or drawCellColor) {
        shader = colorShader;
        shader->setUniform("blending", GLSLValue(mesh->blending()));
  
        glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, (drawCellColor ? cellColorVAid : triangleColorVAid));
        glColorPointer(4, GL_UNSIGNED_BYTE, 0, 0);
      } else {
        setupSurfColorMap();
        setupHeatColorMap();
        setupLabelColorMap();
  
        if(drawImageTex) {
          shader = textureShader;
          glColor4f(0, 0, 0, 1);
          bind2DTex(imgTexId);
    
          shader->setUniform("tex", GLSLValue(Shader::AT_TEX2D));
          shader->setUniform("blending", GLSLValue(mesh->blending()));
    
          glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, imgVAid);
          glfuncs->glTexCoordPointer(2, GL_FLOAT, 0, 0);
        } else if(drawStackTex) {
          shader = volumeShader;
          setup3DRenderingData(shader);
          shader->setUniform("blending", GLSLValue(mesh->blending()));
    
          // Scale texture coordinate!
          glMatrixMode(GL_TEXTURE);
          glPushMatrix();
          glMatrixMode(GL_MODELVIEW);
    
          glColor4f(1, 1, 1, 1);
          glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, posVAid);
          glfuncs->glTexCoordPointer(3, GL_FLOAT, 0, 0);
        } else {
          // Handle signal and label
          Point2f &signalBounds = mesh->signalBounds();
          if(signalBounds[1] <= signalBounds[0])
            signalBounds[1] = signalBounds[0] + 1;
          Point2f &heatMapBounds = mesh->heatMapBounds();
          if(heatMapBounds[1] <= heatMapBounds[0])
            heatMapBounds[1] = heatMapBounds[0] + 1;
					Point2f &triangleValueBounds = mesh->triangleValueBounds();
					if(triangleValueBounds[1] <= triangleValueBounds[0])
						triangleValueBounds[1] = triangleValueBounds[0] + 1;
          // The normal case (including labels and heatmap)
          shader = labelShader;
          shader->setUniform("labels", GLSLValue(mesh->showLabels()));
          shader->setUniform("heatmap", GLSLValue(mesh->showHeat() or mesh->showWallHeat()));
          // shader->setUniform("surfcolormap", GLSLValue(Shader::AT_SURF_TEX));
          // shader->setUniform("labelcolormap", GLSLValue(Shader::AT_LABEL_TEX));
          // shader->setUniform("heatcolormap", GLSLValue(Shader::AT_HEAT_TEX));
          shader->setUniform("blending", GLSLValue(mesh->blending()));
          shader->setUniform("signalBounds", GLSLValue(signalBounds));
          shader->setUniform("heatBounds", GLSLValue(heatMapBounds));
          shader->setUniform("infinity", GLSLValue(std::numeric_limits<float>::infinity()));
    
          glColor4f(1, 1, 1, 1);
          glMatrixMode(GL_TEXTURE);
          glLoadIdentity();
          glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, texVAid);
          glfuncs->glTexCoordPointer(3, GL_FLOAT, 0, 0);
          setupHeatColorMap();
    
          glMatrixMode(GL_MODELVIEW);
        }
      }
    }
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);
  
    if(shader and !select) {
      shader->setUniform("brightness", GLSLValue(powf(mesh->brightness(), 2.0f)));
      shader->setUniform("opacity", GLSLValue(mesh->opacity()));
      if(mesh->surfVertexView() == "Stack Texture") {
        shader->setUniform("second_brightness", GLSLValue(powf(mesh->brightness(), 2.0f)));
        shader->setUniform("second_opacity", GLSLValue(mesh->opacity()));
      }
      const BoundingBox3f& bbox = mesh->boundingBox();
      float meshSize = norm(bbox[1] - bbox[0]);
      shader->setUniform("shift", GLSLValue(meshShift * meshSize));
      shader->useShaders();
      shader->setupUniforms();
      setupSurfColorMap();
    }
    glDrawArrays(GL_TRIANGLES, 0, posVA.size());
  
    if(select)
      glDisableClientState(GL_COLOR_ARRAY);
    else {
      if(shader)
        shader->stopUsingShaders();
      glDisableClientState(GL_NORMAL_ARRAY);
      if(drawTriColor or drawCellColor)
        glDisableClientState(GL_COLOR_ARRAY);
      else
        glDisableClientState(GL_TEXTURE_COORD_ARRAY);
    }
    glDisableClientState(GL_VERTEX_ARRAY);

    // Unbind select framebuffer
    if(select)
      stopSelectFbo();
  
    if(mesh->surfVertexView() == "Stack Texture") {
      glMatrixMode(GL_TEXTURE);
      glPopMatrix();
    }
  
    glMatrixMode(GL_MODELVIEW);
    glPopMatrix();
    unbind3DTex();
    unbind2DTex();
  
    glDisable(GL_POLYGON_OFFSET_FILL);
    glDisable(GL_POLYGON_OFFSET_LINE);
    REPORT_GL_ERROR("drawSurf");
  }

  // Reset mesh
  void ImgData::resetMesh()
  {
    mesh->reset();
    fillVBOs();
  }

  // Reset stack data
  void ImgData::resetStack()
  {
    stack->reset();
    unloadTex();
    emit stackUnloaded();
  }
  
  // Initialize controls
  void ImgData::initControls(QWidget* viewer)
  {
    // Transfer functions
    connect(workTransferDlg, SIGNAL(changedTransferFunction(const TransferFunction &)), this,
            SLOT(updateWorkColorMap(const TransferFunction &)));
    connect(mainTransferDlg, SIGNAL(changedTransferFunction(const TransferFunction &)), this,
            SLOT(updateMainColorMap(const TransferFunction &)));
    connect(surfTransferDlg, SIGNAL(changedTransferFunction(const TransferFunction &)), this,
            SLOT(updateSurfColorMap(const TransferFunction &)));
    connect(heatTransferDlg, SIGNAL(changedTransferFunction(const TransferFunction &)), this,
            SLOT(updateHeatColorMap(const TransferFunction &)));
  
    connect(workTransferDlg, SIGNAL(changedTransferFunction(const TransferFunction &)), viewer, SLOT(UpdateSlot()));
    connect(mainTransferDlg, SIGNAL(changedTransferFunction(const TransferFunction &)), viewer, SLOT(UpdateSlot()));
    connect(surfTransferDlg, SIGNAL(changedTransferFunction(const TransferFunction &)), viewer, SLOT(UpdateSlot()));
    connect(heatTransferDlg, SIGNAL(changedTransferFunction(const TransferFunction &)), viewer, SLOT(UpdateSlot()));
  }
  
  void ImgData::unloadTex()
  {
    if(mainDataTexId)
      glfuncs->glDeleteTextures(1, &mainDataTexId);
    if(workDataTexId)
      glfuncs->glDeleteTextures(1, &workDataTexId);
    if(mcmapTexId)
      glfuncs->glDeleteTextures(1, &mcmapTexId);
    if(wcmapTexId)
      glfuncs->glDeleteTextures(1, &wcmapTexId);
    if(surfTexId)
      glfuncs->glDeleteTextures(1, &surfTexId);
    if(heatTexId)
      glfuncs->glDeleteTextures(1, &heatTexId);
    if(colorMapTexId)
      glDeleteTextures(1, &heatTexId);
  
    surfTexId = heatTexId = colorMapTexId = wcmapTexId = mcmapTexId = mainDataTexId = workDataTexId = 0;
  
    // Histograms require update
    invalidateHistogram(MainHist);
    invalidateHistogram(WorkHist);
    invalidateHistogram(SurfHist);
    invalidateHistogram(HeatHist);
  }
  
  // Get decimation factors for each dimension
  Point3u ImgData::getTexStep()
  {
    Point3u step;
    for(uint i = 0; i < 3; i++) {
      step[i] = stack->size()[i] / MaxTexSize[i];
      if(stack->size()[i] % MaxTexSize[i] > 0)
        step[i]++;
    }
    return (step);
  }
  
  // Load data into texture, handle decimation if texture too big
  void ImgData::load3DTexData(const GLuint texId, const Point3u size, const ushort* data)
  {
    if(swapTextureBytes(texId))
      glPixelStorei(GL_UNPACK_SWAP_BYTES, GL_TRUE);
    // Find if any dimension need to be decimated
    Point3u step = getTexStep();
    if(step == Point3u(1u, 1u, 1u)) {
      GLsizei width = size.x();
  
      if(width % 4 == 0)
        glPixelStorei(GL_UNPACK_ALIGNMENT, 8);
      else if(width % 2 != 0)
        glPixelStorei(GL_UNPACK_ALIGNMENT, 2);
  
      glfuncs->glTexImage3D(GL_TEXTURE_3D, 0, internalFormat(texId), size.x(), size.y(), size.z(), 0, GL_ALPHA,
			    GL_UNSIGNED_SHORT, data);
    } else {
      Point3u sz = Point3u(size.x() / step.x(), size.y() / step.y(), size.z() / step.z());
  
      HVecUS tpix(size_t(sz.x()) * sz.y() * sz.z());
      ushort* p = tpix.data();
  
      // Decimate array
      for(uint z = step.z() / 2; z < sz.z() * step.z(); z += step.z())
        for(uint y = step.y() / 2; y < sz.y() * step.y(); y += step.y())
          for(uint x = step.x() / 2; x < sz.x() * step.x(); x += step.x())
            *p++ = data[offset(x, y, z)];
  
      GLsizei width = sz.x();
  
      if(width % 4 == 0)
        glPixelStorei(GL_UNPACK_ALIGNMENT, 8);
      else if(width % 2 != 0)
        glPixelStorei(GL_UNPACK_ALIGNMENT, 2);
  
      // Load smaller texture
      glfuncs->glTexImage3D(GL_TEXTURE_3D, 0, internalFormat(texId), sz.x(), sz.y(), sz.z(), 0, GL_ALPHA, GL_UNSIGNED_SHORT,
			    tpix.data());
    }
    glPixelStorei(GL_UNPACK_ALIGNMENT, 4);   // Restore default
    if(swapTextureBytes(texId))
      glPixelStorei(GL_UNPACK_SWAP_BYTES, GL_FALSE);
  
    // Histograms require update
    if(texId == mainDataTexId)
      invalidateHistogram(MainHist);
    else if(texId == workDataTexId)
      invalidateHistogram(WorkHist);
  
    glfuncs->glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, interpolation(mainDataTexId));
    glfuncs->glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, interpolation(mainDataTexId));
    glfuncs->glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_BORDER);
    glfuncs->glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_BORDER);
    glfuncs->glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_BORDER);
    float anisotropy;
    glfuncs->glGetFloatv(GL_MAX_TEXTURE_MAX_ANISOTROPY_EXT, &anisotropy);
    glfuncs->glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_MAX_ANISOTROPY_EXT, anisotropy);
  }
  
  void ImgData::reloadTex(GLuint texId)
  {
    if(texId == 0 or stack->storeSize() == 0
       or (texId == mainDataTexId and stack->main()->data().size() != stack->storeSize())
       or (texId == workDataTexId and stack->work()->data().size() != stack->storeSize())) {
      unloadTex();
      return;
    }
  
    glBindTexture(GL_TEXTURE_3D, texId);
    HVecUS& wdata = stack->work()->data();
    HVecUS& mdata = stack->main()->data();
    const ushort* data = (texId == mainDataTexId) ? mdata.data() : wdata.data();
    load3DTexData(texId, stack->size(), data);
  }
  
  // Update a section of a texture defined in image coordinates by a bounding box
  void ImgData::updateTex(GLuint texId, BoundingBox3i bBox)
  {
    if(texId == 0 or stack->storeSize() == 0 or stack->main()->data().size() != stack->storeSize()
       or stack->work()->data().size() != stack->storeSize())
      return;
  
    HVecUS& data = (texId == mainDataTexId ? stack->main()->data() : stack->work()->data());
  
    Point3i base(bBox[0].x(), bBox[0].y(), bBox[0].z());
    Point3i bsz(bBox[1].x() - bBox[0].x(), bBox[1].y() - bBox[0].y(), bBox[1].z() - bBox[0].z());
  
    if(swapTextureBytes(texId))
      glPixelStorei(GL_UNPACK_SWAP_BYTES, GL_TRUE);
    // Get factor to decimate texture if required
    Point3u step = getTexStep();
    // Enlarge bounding box to multiples of steps
    if(step != Point3u(1u, 1u, 1u)) {
      base.x() -= base.x() % step.x();
      base.y() -= base.y() % step.y();
      base.z() -= base.z() % step.z();
      bsz.x() += base.x() % step.x();
      bsz.y() += base.y() % step.y();
      bsz.z() += base.z() % step.z();
      if(bsz.x() % step.x() > 0)
        bsz.x() = bsz.x() / step.x() + 1;
      else
        bsz.x() /= step.x();
      if(bsz.y() % step.y() > 0)
        bsz.y() = bsz.y() / step.y() + 1;
      else
        bsz.y() /= step.y();
      if(bsz.z() % step.z() > 0)
        bsz.z() = bsz.z() / step.z() + 1;
      else
        bsz.z() /= step.z();
  
      HVecUS tpix(size_t(bsz.x()) * bsz.y() * bsz.z());
  
      ushort* p = tpix.data();
      for(uint z = base.z() + step.z() / 2; z < base.z() + bsz.z() * step.z(); z += step.z())
        for(uint y = base.y() + step.y() / 2; y < base.y() + bsz.y() * step.y(); y += step.y())
          for(uint x = base.x() + step.x() / 2; x < base.x() + bsz.x() * step.x(); x += step.x())
            *p++ = data[offset(x, y, z)];
  
      glPixelStorei(GL_UNPACK_ALIGNMENT, 2);
  
      // Load texture subimage
      glfuncs->glBindTexture(GL_TEXTURE_3D, texId);
      glfuncs->glTexSubImage3D(GL_TEXTURE_3D, 0, base.x() / step.x(), base.y() / step.y(), base.z() / step.z(), bsz.x(),
			       bsz.y(), bsz.z(), GL_ALPHA, GL_UNSIGNED_SHORT, tpix.data());
      glfuncs->glBindTexture(GL_TEXTURE_3D, 0);
  
      glfuncs->glPixelStorei(GL_UNPACK_ALIGNMENT, 4);     // Restore default
    } else {
      // Update only the part needed
      glfuncs->glPixelStorei(GL_UNPACK_ALIGNMENT, 2);
  
      glfuncs->glBindTexture(GL_TEXTURE_3D, texId);
      glfuncs->glPixelStorei(GL_UNPACK_ROW_LENGTH, stack->size().x());
      glfuncs->glPixelStorei(GL_UNPACK_IMAGE_HEIGHT, stack->size().y());
      glfuncs->glTexSubImage3D(GL_TEXTURE_3D, 0, base.x(), base.y(), base.z(), bsz.x(), bsz.y(), bsz.z(), GL_ALPHA,
			       GL_UNSIGNED_SHORT, &data[offset(base.x(), base.y(), base.z())]);
      glfuncs->glPixelStorei(GL_UNPACK_ROW_LENGTH, 0);
      glfuncs->glPixelStorei(GL_UNPACK_IMAGE_HEIGHT, 0);
      glfuncs->glBindTexture(GL_TEXTURE_3D, 0);
  
      glPixelStorei(GL_UNPACK_ALIGNMENT, 4);     // Restore default
    }
    if(swapTextureBytes(texId))
      glPixelStorei(GL_UNPACK_SWAP_BYTES, GL_FALSE);
  
    // Histograms require update
    if(texId == mainDataTexId)
      invalidateHistogram(MainHist);
    else if(texId == workDataTexId)
      invalidateHistogram(WorkHist);

    REPORT_GL_ERROR("updateTex");
  }
  
  void ImgData::reloadMainTex(const BoundingBox3i& bbox)
  {
    if(mainDataTexId) {
      reloadTex(mainDataTexId);
      if(bbox.empty() or bbox.size() == Point3i(stack->size())) {
        reloadTex(mainDataTexId);
      } else {
        updateTex(mainDataTexId, bbox);
      }
    } else {
      if(stack->storeSize() == 0 or stack->main()->data().size() != stack->storeSize()) {
        Information::out << "Error cannot reload main tex as SizeXYZ = " << stack->storeSize()
                         << " and MainData.size() = " << stack->main()->data().size() << endl;
        unloadTex();
        return;
      }
      Information::out << "Creating new texture of size " << stack->size() << endl;
      glGenTextures(1, &mainDataTexId);
      glBindTexture(GL_TEXTURE_3D, mainDataTexId);
      HVecUS& mdata = stack->main()->data();
      load3DTexData(mainDataTexId, stack->size(), mdata.data());
    }
  }
  
  void ImgData::reloadWorkTex(const BoundingBox3i& bbox)
  {
    if(workDataTexId) {
      if(bbox.empty() or bbox.size() == Point3i(stack->size())) {
        reloadTex(workDataTexId);
      } else {
        updateTex(workDataTexId, bbox);
      }
    } else {
      if(stack->storeSize() == 0 or stack->work()->data().size() != stack->storeSize()) {
        Information::out << "Error cannot reload work tex as SizeXYZ = " << stack->storeSize()
                         << " and WorkData.size() = " << stack->main()->data().size() << endl;
        unloadTex();
        return;
      }
      glGenTextures(1, &workDataTexId);
      glBindTexture(GL_TEXTURE_3D, workDataTexId);
      HVecUS& wdata = stack->work()->data();
      load3DTexData(workDataTexId, stack->size(), wdata.data());
    }
  }
  
  void ImgData::loadTex()
  {
    if(stack->storeSize() == 0 or stack->main()->data().size() != stack->storeSize()
       or stack->work()->data().size() != stack->storeSize()) {
      unloadTex();
      return;
    }
  
    // Main stack
    if(!mainDataTexId)
      glGenTextures(1, &mainDataTexId);
    glBindTexture(GL_TEXTURE_3D, mainDataTexId);
    HVecUS& mdata = stack->main()->data();
    load3DTexData(mainDataTexId, stack->size(), mdata.data());
    // Work stack
    if(!workDataTexId)
      glGenTextures(1, &workDataTexId);
    glBindTexture(GL_TEXTURE_3D, workDataTexId);
    HVecUS& wdata = stack->work()->data();
    load3DTexData(workDataTexId, stack->size(), wdata.data());
  }
  
  // Bind the stack texture
  void ImgData::bind3DTex(GLuint texId, Shader::ActiveTextures atexId)
  {
    if(texId) {
      glMatrixMode(GL_TEXTURE);
      glLoadIdentity();
      const Point3f& Size = multiply(stack->step(), Point3f(stack->size()));
      glScaled(1 / Size.x(), 1 / Size.y(), 1 / Size.z());
      Point3f Origin = stack->origin();
      glTranslatef(-Origin.x(), -Origin.y(), -Origin.z());
      Shader::activeTexture(atexId);
      glBindTexture(GL_TEXTURE_3D, texId);
      glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, interpolation(texId));
      glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, interpolation(texId));
      Shader::activeTexture(Shader::AT_NONE);
      glBindTexture(GL_TEXTURE_3D, 0);
      glMatrixMode(GL_MODELVIEW);
    }
  }
  
  // Bind the stack texture
  void ImgData::bind2DTex(GLuint texId)
  {
    if(texId) {
      glMatrixMode(GL_TEXTURE);
      glLoadIdentity();
      Shader::activeTexture(Shader::AT_TEX2D);
      glBindTexture(GL_TEXTURE_2D, texId);
      glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, interpolation(texId));
      glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, interpolation(texId));
      Shader::activeTexture(Shader::AT_NONE);
      glBindTexture(GL_TEXTURE_2D, 0);
      glMatrixMode(GL_MODELVIEW);
    }
  }
  
  void ImgData::unbind3DTex()
  {
    glMatrixMode(GL_TEXTURE);
    glLoadIdentity();
    Shader::activeTexture(Shader::AT_TEX3D);
    glBindTexture(GL_TEXTURE_3D, 0);
    Shader::activeTexture(Shader::AT_SECOND_TEX3D);
    glBindTexture(GL_TEXTURE_3D, 0);
    Shader::activeTexture(Shader::AT_NONE);
    glMatrixMode(GL_MODELVIEW);
  }
  
  void ImgData::unbind2DTex()
  {
    glMatrixMode(GL_TEXTURE);
    glLoadIdentity();
    Shader::activeTexture(Shader::AT_TEX2D);
    glBindTexture(GL_TEXTURE_2D, 0);
    Shader::activeTexture(Shader::AT_NONE);
    glMatrixMode(GL_MODELVIEW);
  }
  
  // Load a texture image (for Keyence)
  void ImgData::loadImgTex(const QImage& image)
  {
    if(image.isNull() or image.width() <= 0 or image.height() <= 0)
      return;
    if(imgTexId)
      glfuncs->glDeleteTextures(1, &imgTexId);
    glGenTextures(1, &imgTexId);
    glBindTexture(GL_TEXTURE_2D, imgTexId);
  
    std::vector<float> pix(image.width() * image.height() * 3);
    float* px = &pix[0];
    for(int y = image.height() - 1; y >= 0; y--)
      for(int x = 0; x < image.width(); x++) {
        *px++ = float(qRed(image.pixel(x, y))) / 255.0;
        *px++ = float(qGreen(image.pixel(x, y))) / 255.0;
        *px++ = float(qBlue(image.pixel(x, y))) / 255.0;
      }
  
    glfuncs->glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, image.width(), image.height(), 0,
			  GL_RGB, GL_FLOAT, &pix[0]);
  
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_BORDER);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_BORDER);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_BORDER);
  }
  
  void ImgData::initTex()
  {
    // Initialize 1D index texture
    LabelColorsChanged = true;
    reloadLabelTex();
  }
  
  void ImgData::reloadLabelTex()
  {
    if(LabelColorsChanged) {
      if(!labelTexId)
        glGenTextures(1, &labelTexId);
      glBindTexture(GL_TEXTURE_1D, labelTexId);
  
      if(!LabelColors.empty())
        glfuncs->glTexImage1D(GL_TEXTURE_1D, 0, GL_RGB, LabelColors.size(), 0,
			      GL_RGBA, GL_FLOAT, &LabelColors[0]);
      else
        glfuncs->glTexImage1D(GL_TEXTURE_1D, 0, GL_RGB, 0, 0, GL_RGB, GL_FLOAT, 0);
      glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
      glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
      glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_REPEAT);
      glBindTexture(GL_TEXTURE_1D, 0);
      LabelColorsChanged = false;
    }
  }
  
  void ImgData::correctSelection(bool inclusive) 
  {
    mesh->correctSelection(inclusive);
  }
  
  // Select vertices with a specified label
  // Repeat used for "select on mesh" option on "Edit Labels" dialog
  void ImgData::selectLabel(int label, int repeat)
  {
    vvGraph& S = mesh->graph();
    if(!mesh->showSurface() or S.empty())
      return;
  
    bool changed = false;
    if(repeat > 0) {
      #pragma omp parallel for
      for(uint i = 0; i < S.size(); i++) {
        vertex v = S[i];
        if((v->label % repeat) == label and !v->selected) {
          changed = true;
          v->selected = true;
        }
      }
    } else {
      #pragma omp parallel for
      for(uint i = 0; i < S.size(); i++) {
        vertex v = S[i];
        if(v->label == label and !v->selected) {
          changed = true;
          v->selected = true;
        }
      }
    }
    if(changed) {
      mesh->correctSelection(true);
      updateSelection();
    }
  }

  // Select vertices from a set of labels
  void ImgData::selectLabel(const IntSet &labels)
  {
    vvGraph& S = mesh->graph();
    if(labels.empty() or S.empty())
      return;
  
    bool changed = false;
    #pragma omp parallel for
    for(uint i = 0; i < S.size(); i++) {
      vertex v = S[i];
      if(labels.count(v->label) > 0 and !v->selected) {
        changed = true;
        v->selected = true;
      }
    }
    if(changed) {
      mesh->correctSelection(true);
      updateSelection();
    }
  } 

  // Unselect vertices by label
  void ImgData::unselectLabel(const IntSet &labels)
  {
    vvGraph& S = mesh->graph();
    if(labels.empty() or S.empty())
      return;
  
    bool changed = false;
    #pragma omp parallel for
    for(uint i = 0; i < S.size(); i++) {
      vertex v = S[i];
      if(labels.count(v->label) > 0 and v->selected) {
        changed = true;
        v->selected = false;
      }
    }
    if(changed) {
      mesh->correctSelection(true);
      updateSelection();
    }
  }
  
  // Select vertices by parent
  void ImgData::selectParent(int parent, int repeat)
  {
    vvGraph& S = mesh->graph();
    if(!mesh->showSurface() or S.empty())
      return;
  
    bool changed = false;
    if(repeat > 0) {
      const IntIntAttr& parentMap = mesh->parents();
      #pragma omp parallel for
      for(uint i = 0; i < S.size(); i++) {
        vertex v = S[i];
        if(v->selected)
          continue;
        IntIntAttr::const_iterator it = parentMap.find(v->label);
        if(it != parentMap.end() and it->second % repeat == parent) {
          changed = true;
          v->selected = true;
        }
      }
    } else {
      const IntIntAttr& parentMap = mesh->parents();
      #pragma omp parallel for
      for(uint i = 0; i < S.size(); i++) {
        vertex v = S[i];
        if(v->selected)
          continue;
        IntIntAttr::const_iterator it = parentMap.find(v->label);
        if(it != parentMap.end() and it->second == parent) {
          changed = true;
          v->selected = true;
        }
      }
    }
    if(changed) {
      mesh->correctSelection(true);
      updateSelection();
    }
  }

  // Select vertices by parent
  void ImgData::selectParent(const IntSet &parents)
  {
    vvGraph& S = mesh->graph();
    if(parents.empty() or S.empty())
      return;
  
    bool changed = false;
    const IntIntAttr& parentMap = mesh->parents();
    #pragma omp parallel for
    for(uint i = 0; i < S.size(); i++) {
      vertex v = S[i];
      if(v->selected)
        continue;
      IntIntAttr::const_iterator it = parentMap.find(v->label);
      if(it != parentMap.end() and parents.count(it->second) > 0) {
        changed = true;
        v->selected = true;
      }
    }
    if(changed) {
      mesh->correctSelection(true);
      updateSelection();
    }
  } 

  // Unselect vertices by
  void ImgData::unselectParent(const IntSet &parents)
  {
    vvGraph& S = mesh->graph();
    if(parents.empty() or S.empty())
      return;
  
    bool changed = false;
    const IntIntAttr& parentMap = mesh->parents();
    #pragma omp parallel for
    for(uint i = 0; i < S.size(); i++) {
      vertex v = S[i];
      if(!v->selected)
        continue;
      IntIntAttr::const_iterator it = parentMap.find(v->label);
      if(it != parentMap.end() and parents.count(it->second) > 0) {
        changed = true;
        v->selected = false;
      }
    }
    if(changed) {
      mesh->correctSelection(true);
      updateSelection();
    }
  }
  
  // Select connected vertices
  void ImgData::selectConnected(std::vector<uint>& vList, bool unselect)
  {
    vvGraph& S = mesh->graph();
    if(!mesh->showSurface() or S.empty() or vList.empty())
      return;
  
    vvGraph V;
    for(uint i = 0; i < vList.size(); i++) {
      vertex v = idVA[vList[i]];
      V.insert(v);
    }
    mesh->getConnectedVertices(S, V);
    bool changed = false;
    #pragma omp parallel for
    for(uint i = 0; i < V.size(); i++) {
      vertex v = V[i];
      if( v->selected != !unselect) {
        changed = true;
        v->selected = !unselect;
      }
    }
    if(changed)
      updateSelection();
  }

  // Select connected vertices
  void ImgData::selectVertices(std::vector<uint>& vList, bool unselect)
  {
    vvGraph& S = mesh->graph();
    if(!mesh->showSurface() or S.empty() or vList.empty())
      return;
  
    bool changed = false;
    for(uint i = 0; i < vList.size(); i++) {
      vertex v = idVA[vList[i]];
      if(v->selected != !unselect) {
        v->selected = !unselect;
        changed = true;
      }
    }
    if(changed)
      updateSelection();
  } 

  void ImgData::fillLabel(const IntSet &labels, int currlabel)
  {
    vvGraph& S = mesh->graph();
    if(labels.empty() or S.empty())
      return;

    mesh->setLabel(std::max(currlabel, mesh->viewLabel()));

    #pragma omp parallel for
    for(uint i = 0; i < S.size(); i++) {
      vertex v = S[i];
      if(labels.count(v->label) > 0)
        v->label = currlabel;
    }
  
    // Clear cell axis, cell center for new label
    forall(int label, labels) {
      mesh->labelCenter().erase(label);
      mesh->labelNormal().erase(label);
      mesh->labelCenterVis().erase(label);
      mesh->labelNormalVis().erase(label);
      mesh->cellAxis().erase(label);
      mesh->cellAxisVis().erase(label);
      mesh->cellAxisColor().erase(label);
      // Clear cell axis, cell center for current label
      mesh->labelCenter().erase(currlabel);
      mesh->labelNormal().erase(currlabel);
      mesh->labelCenterVis().erase(currlabel);
      mesh->labelNormalVis().erase(currlabel);
      mesh->cellAxis().erase(currlabel);
      mesh->cellAxisVis().erase(currlabel);
      mesh->cellAxisColor().erase(currlabel);
      // Clear parent label
      IntIntAttr::iterator p = mesh->parents().find(label);
      if(p != mesh->parents().end()) {
        int lparent = p->second;
        mesh->parentCenter().erase(lparent);
        mesh->parentNormal().erase(lparent);
        mesh->parentCenterVis().erase(lparent);
        mesh->parentNormalVis().erase(lparent);
        mesh->parents().erase(label);
      }
    }
    // Clear parent label
    IntIntAttr::iterator cp = mesh->parents().find(currlabel);
    if(cp != mesh->parents().end()) {
      int currparent = cp->second;
      mesh->parentCenter().erase(currparent);
      mesh->parentNormal().erase(currparent);
      mesh->parentCenterVis().erase(currparent);
      mesh->parentNormalVis().erase(currparent);
      mesh->parents().erase(currlabel);
    }
  
    // Update margin
    markMargins(S);
    updateTriColor();
    updateLines();
  }

  void ImgData::setParent(int label, int parentLabel)
  {
    IntSet labels;
    labels.insert(label);
    setParent(labels, parentLabel);
  }

  void ImgData::setParent(const IntSet &labels, int parentLabel)
  {
    bool changed = false;
  
    IntIntAttr & parents = mesh->parents();
    forall(int label, labels) {
      if(parentLabel == 0) {
        parents.erase(label);
        changed = true;
      } else if(parents[label] != parentLabel) {
        parents[label] = parentLabel;
        changed = true;
      }
    }
  
    if(changed)
      updateTriColor();
  } 

  void ImgData::addSeed(int label, const std::vector<uint> &vList)
  {
    if(vList.empty())
      return;
  
    mesh->setLabel(std::max(label, mesh->viewLabel()));

    if(mesh->useParents()) { // if parents are activated, assign parent labels to connected regions (3d cells)
      vvGraph& S = mesh->graph();
      vvGraph V;
      for(uint i = 0; i < vList.size(); i++) {
        vertex v = idVA[vList[i]];
        V.insert(v);
      }
      mesh->getConnectedVertices(S, V);

      // prevent overwriting all parents in a 2.5D mesh
      if(S.size() == V.size()){
        Information::out << "Can't add parents in a non-volumetric mesh. Deactivate parents or use a mesh with volumetric cells" << endl;
        return;
      }

      #pragma omp parallel for
      for(uint i = 0; i < V.size(); i++) {
        vertex v = V[i];
        mesh->parents()[v->label] = label;
      }
      updateTriColor();
    } else { // if not, add seeds
      // Label by cell
      if(mesh->meshType() == "MGXC") {
        std::set<uint> vset;
        forall(uint u, vList) {
          vertex v = idVA[u];
          if(v->type == 'c') {
            vset.insert(u);

            // Find all triangles that contain the center vertex
            for(uint i = 0; i < selVA.size(); i += 3)
              for(uint j = 0; j < 3; j++)
                if(idVA[i + j] == idVA[u])
                  for(uint k = 0; k < 3; k++)
                    vset.insert(i + k);
          }
        }
        std::vector<uint> tlist;
        forall(uint u, vset)
          tlist.push_back(u);
        updLabel(label, tlist);
      } else if(mesh->meshType() == "MGXM")
        updLabel(label, vList);
    }
  }
 
  void ImgData::drawSignal(int signal, const std::vector<uint> &vList)
  {
    updSignal(signal, vList);
  }
   
  // Fill selection with label
  void ImgData::fillSelect(int label)
  {
    vvGraph& S = mesh->graph();
    if(S.empty() or !mesh->showSurface())
      return;
  
    if(mesh->useParents()) {
      IntIntAttr &parents = mesh->parents();
      #pragma omp parallel for
      for(uint i = 0; i < selectV.size(); i++) {
        vertex v = selectV[i];
        parents[v->label] = label;
      }
    } else {
      #pragma omp parallel for
      for(uint i = 0; i < selectV.size(); i++) {
        vertex v = selectV[i];
        v->label = label;
      }
		}
  
    markMargins(S);
    updateTriColor();
    updateLines();
  }
  
  // Delete Label
  void ImgData::deleteLabel(int label)
  {
    vvGraph& S = mesh->graph();
    if(S.empty())
      return;
  
    std::vector<vertex> T;
  
    if(MeshSelect) {
      if(mesh->showMeshPoints()) {
        forall(const vertex &v, selectV)
          T.push_back(v);
        selectV.clear();
      }
    } else if(mesh->showSurface() and mesh->showLabels()) {
      // Label delete, remove all vertices with this label
      // or all border vertices connected only to other border vertices or this label
      forall(const vertex& v, S) {
        if(v->label == label)
          T.push_back(v);
        else if(v->label == -1) {
          bool remove = true;
          forall(const vertex& n, S.neighbors(v)) {
            if(n->label != label and n->label != -1) {
              remove = false;
              break;
            }
          }
          if(remove)
            T.push_back(v);
        }
      }
    } else
      return;
  
    // Now delete them
    forall(const vertex& v, T)
      S.erase(v);
  
    // Delete any singletons
    T.clear();
    forall(const vertex& v, S)
      if(S.valence(v) <= 1)
        T.push_back(v);
    forall(const vertex& v, T)
      S.erase(v);
  
    fillVBOs();
    if(MeshSelect and mesh->showMeshPoints())
      Information::setStatus(QString("%1 deleted selected vertices").arg(Section));
    else
      Information::setStatus(QString("%1 deleted label:").arg(Section).arg(label));
  }

  // Update only the position of points and lines
  void ImgData::updatePos()
  { 
    // Bounding box
    BoundingBox3f bBox(Point3f(FLT_MAX, FLT_MAX, FLT_MAX), Point3f(-FLT_MAX, -FLT_MAX, -FLT_MAX));

    vvGraph &S = mesh->graph();
    if(pntsVA.size() != S.size()) {
      Information::out << "updatePos:Error, point count mismatch in pntds vertex array" << endl;
      return;
    }   
    if(S.empty())
      return;

    for(size_t i = 0; i < S.size(); i++)
      bBox |= pntsVA[i] = Point3f(S[i]->pos);
    mesh->setBoundingBox(bBox);
 
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, pntsVAid);
    glfuncs->glBufferData(GL_ARRAY_BUFFER_ARB, pntsVA.size() * sizeof(Point3f), &pntsVA[0], GL_STATIC_DRAW_ARB); 
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);

    REPORT_GL_ERROR("updatePos");
  }
 
  // Update points for lines and calculate vertex, line and triangle counts, set saveIds
  void ImgData::updatePoints(uint &VCount, uint &LCount, uint &TCount)
  { 
    uint badEdges = 0;
    VCount = LCount = TCount = 0;

    // Clear vertex arrays
    idVA.clear();
    triCell.clear();

    // Bounding box
    BoundingBox3f bBox(Point3f(FLT_MAX, FLT_MAX, FLT_MAX), Point3f(-FLT_MAX, -FLT_MAX, -FLT_MAX));

    vvGraph &S = mesh->graph();
    pntsVA.resize(S.size());
    if(S.empty())
      return;

    bool mesh3D = mesh->mesh3D();

    // We use S graph for lines for both 2D and 3D
    forall(const vertex &v, S) {
      pntsVA[VCount] = Point3f(v->pos);
      v->saveId = VCount++;

      bBox |= Point3f(v->pos);

      forall(const vertex& n, S.neighbors(v)) {
        if(!S.edge(n, v)) {
          badEdges++;
          n->selected = v->selected = true;
        }
        if(S.uniqueLine(v, n))
          LCount++;
        vertex m = S.nextTo(v, n);

        if(!mesh3D and S.uniqueTri(v, n, m)) {
          idVA.push_back(v); 
          idVA.push_back(n); 
          idVA.push_back(m); 
        }
      }
    }
    // Faces need to use 3D graph
    if(mesh3D) {
      forall(const cell &c, mesh->cells())
        forall(const vertex &v, c->S)
          forall(const vertex& n, c->S.neighbors(v)) {
            vertex m = c->S.nextTo(v, n);

            if(c->S.uniqueTri(v, n, m)) {
              idVA.push_back(v); 
              idVA.push_back(n); 
              idVA.push_back(m); 
              triCell.push_back(c);
            }
          }
    }
    TCount = idVA.size()/3;
    mesh->setBoundingBox(bBox);
 
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, pntsVAid);
    glfuncs->glBufferData(GL_ARRAY_BUFFER_ARB, pntsVA.size() * sizeof(Point3f), &pntsVA[0], 
                                                                       GL_STATIC_DRAW_ARB); 
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);

    if(badEdges)
      Information::out << badEdges << " bad edges found in graph." << endl;

    REPORT_GL_ERROR("updatePoints");
  }
 
  void ImgData::updateNextLabel()
  {
    int nextLabel = mesh->viewLabel();
    vvGraph &S = mesh->graph();
    forall(const vertex &v, S)
      if(v->label > nextLabel)
        nextLabel = v->label; 
    if(nextLabel > mesh->viewLabel())
      mesh->setLabel(nextLabel);
  }
      
  void ImgData::updateTriColor()
  {
    if(idVA.empty())
      return;

    updateNextLabel();

    if(idVA.size() != texVA.size()) {
      Information::out << "updateTriColor:Error, point count mismatch in tex vertex array" << endl;
      return;
    }   

    bool mesh3D = mesh->mesh3D();
    if(mesh3D and idVA.size()/3 != triCell.size()) {
      Information::out << "updateTriColor:Error, triCell size incorrect" << endl;
      return;
    }   
		uint triangles = idVA.size()/3;
    const IntIntAttr &parents = mesh->parents();
    const IntFloatAttr &labelHeat = mesh->labelHeat();
    bool drawWallHeat = mesh->surfLabelView() == "Wall Heat";
    IntIntFloatAttr &wallHeat = mesh->wallHeat();
    
    TriFloatAttr &triangleValue = mesh->triangleValue();
    TriIntAttr &triangleIndex = mesh->triangleIndex();
    TriVec3ColorbAttr &triangleColor = mesh->triangleColor();

    ColorMap triValueColorMap(&mesh->triangleValueColors(), &mesh->triangleValueBounds());
    ColorMap triIndexColorMap(&mesh->triangleIndexColors());

    bool drawTriValue = mesh->surfTriangleView() == "Triangle Value" and triangleValue.size() > 0;
    bool drawTriIndex = mesh->surfTriangleView() == "Triangle Index" and triangleIndex.size() > 0;
    bool drawTriColor = mesh->surfTriangleView() == "Triangle Color" and triangleColor.size() > 0;
    bool drawTri = drawTriValue or drawTriIndex or drawTriColor;
    triangleColorVA.resize(drawTri ? idVA.size() : 0);

    CellTissue &T = mesh->tissue();
    bool drawCellColor = mesh->surfLabelView() == "Cell Color" and T.hasCellGraph();
    cellColorVA.resize(drawCellColor ? idVA.size() : 0);
    CellVec2ColorbAttr &cellColor = mesh->cellColor();

    IntCellMap labelToCell;
    if(drawCellColor) {
      forall(const cell &c, T.C)
        labelToCell[c->label] = c;
    }
    // Put signal and heat in VBOs
    #pragma omp parallel for
    for(uint i = 0; i < triangles; ++i) {
      uint tcnt = i * 3;
      vertex v(idVA[tcnt]), n(idVA[tcnt + 1]), m(idVA[tcnt + 2]);
      // For 3D case need to get this from cell graph
      int label = 0;
      if(mesh3D)
        label = mesh->getLabel(triCell[i]->label, parents);
      else
        label = mesh->getLabel(v, n, m, parents);
      if(drawWallHeat) {
        texVA[tcnt] = texCoordWalls(label, v, v, n, m, wallHeat);
        texVA[tcnt+1] = texCoordWalls(label, n, v, n, m, wallHeat);
        texVA[tcnt+2] = texCoordWalls(label, m, v, n, m, wallHeat);
      } else {
        texVA[tcnt] = texCoordLabels(label, v, labelHeat);
        texVA[tcnt+1] = texCoordLabels(label, n, labelHeat);
        texVA[tcnt+2] = texCoordLabels(label, m, labelHeat);
      }

      if(drawTriValue) {
        	//const Colorb color = triValueColorMap.getColor(triangleValue[Triangle(v,n,m)]);
				Colorb color;
				TriFloatAttr::iterator it = triangleValue.find(Triangle(v,n,m)); 
				if(it != triangleValue.end()) 
        	color = triValueColorMap.getColor(it->second);
				else 
					color = Colorb(127,127,127,255);
        triangleColorVA[tcnt] = color;
        triangleColorVA[tcnt+1] = color;
        triangleColorVA[tcnt+2] = color;
      }
      if(drawTriIndex) {
        const Colorb &color = triIndexColorMap.getColor(triangleIndex[Triangle(v,n,m)]);
        triangleColorVA[tcnt] = color;
        triangleColorVA[tcnt+1] = color;
        triangleColorVA[tcnt+2] = color;
      }
      if(drawTriColor) {
        Vec3Colorb &colors = triangleColor[Triangle(v,n,m)];
        triangleColorVA[tcnt] = colors[0];
        triangleColorVA[tcnt+1] = colors[1];
        triangleColorVA[tcnt+2] = colors[2];
      }
      if(drawCellColor) {
        int label = getLabel(v, n, m);
        if(labelToCell.count(label) == 1) {
          Vec2Colorb colors = cellColor[labelToCell[label]];
          cellColorVA[tcnt] = colors[(v->label > 0 ? 0 : 1)];
          cellColorVA[tcnt+1] = colors[(n->label > 0 ? 0 : 1)];
          cellColorVA[tcnt+2] = colors[(m->label > 0 ? 0 : 1)];
        }
      }
    }

		// Write to GPU
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, texVAid);
    glfuncs->glBufferData(GL_ARRAY_BUFFER_ARB, texVA.size() 
                                       * sizeof(Point3f), &texVA[0], GL_STATIC_DRAW_ARB);
    if(drawTriValue || drawTriIndex || drawTriColor) {
      glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, triangleColorVAid);
      glfuncs->glBufferData(GL_ARRAY_BUFFER_ARB, triangleColorVA.size() 
                                       * sizeof(Colorb), &triangleColorVA[0], GL_STATIC_DRAW_ARB);
    }
    if(drawCellColor) {
      glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, cellColorVAid);
      glfuncs->glBufferData(GL_ARRAY_BUFFER_ARB, cellColorVA.size() 
                                       * sizeof(Colorb), &cellColorVA[0], GL_STATIC_DRAW_ARB);
    }
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);

    REPORT_GL_ERROR("updateTriColor");
  }
  
  void ImgData::updateTriPos()
  {
    if(idVA.size() != posVA.size()) {
      Information::out << "updateTriPos:Error: point count mismatch in pos vertex array" << endl;
      return;
    }  
    if(idVA.size() != nrmlVA.size()) {
      Information::out << "updateTriPos:Error: point count mismatch in nrml vertex array" << endl;
      return;
    } 

    // Update positions for triangles
    uint triangles = idVA.size()/3;
    float step = std::max(stack->step().x(), std::max(stack->step().y(), stack->step().z()));
    if(step <= 0) 
      step = 1.0;
     
    bool mesh3D = mesh->mesh3D();

    #pragma omp parallel for
    for(uint i = 0; i < triangles; ++i) {
      uint tcnt = i * 3;
      vertex v(idVA[tcnt]), n(idVA[tcnt + 1]), m(idVA[tcnt + 2]);

      // For 3D meshes we need to calculate the normals on the fly
      if(mesh3D) {
        Point3d nrml(0,0,0);
        vvGraph &S = triCell[i]->S;

        calcNormal(S, v, nrml);
        posVA[tcnt] = Point3f(v->pos) - SurfOffset * step * Point3f(nrml);
        nrmlVA[tcnt++] = Point3f(nrml);
  
        calcNormal(S, n, nrml);
        posVA[tcnt] = Point3f(n->pos) - SurfOffset * step * Point3f(nrml);
        nrmlVA[tcnt++] = Point3f(nrml);
  
        calcNormal(S, m, nrml);
        posVA[tcnt] = Point3f(m->pos) - SurfOffset * step * Point3f(nrml);
        nrmlVA[tcnt] = Point3f(nrml);
      } else {
        posVA[tcnt] = Point3f(v->pos) - SurfOffset * step * Point3f(v->nrml);
        nrmlVA[tcnt++] = Point3f(v->nrml);
    
        posVA[tcnt] = Point3f(n->pos) - SurfOffset * step * Point3f(n->nrml);
        nrmlVA[tcnt++] = Point3f(n->nrml);
  
        posVA[tcnt] = Point3f(m->pos) - SurfOffset * step * Point3f(m->nrml);
        nrmlVA[tcnt] = Point3f(m->nrml);
      }
    }
  
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, posVAid);
    glfuncs->glBufferData(GL_ARRAY_BUFFER_ARB, posVA.size() * sizeof(Point3f), &posVA[0], 
                                                                          GL_STATIC_DRAW_ARB);
  
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, nrmlVAid);
    glfuncs->glBufferData(GL_ARRAY_BUFFER_ARB, nrmlVA.size() * sizeof(Point3f), &nrmlVA[0], 
                                                                          GL_STATIC_DRAW_ARB);
  
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);

    REPORT_GL_ERROR("updateTriPos");
  }

  // Update point and lines VAs for cells and borders, also updates selection
  void ImgData::updateLines()
  {
    pbrdVA.clear();
    pcellVA.clear();
    pselVA.clear();

    lbrdVA.clear();
    lcellVA.clear();
    lselVA.clear();

    selectV.clear();
    uint lcnt = 0; // line counter

    const vvGraph &S = mesh->graph();
    pcolVA.resize(S.size());
    if(S.empty())
      return;
  
    for(size_t i = 0; i < S.size(); i++) {
      vertex v = S[i];
      // Point/line color
      pcolVA[i] = pColor(v);
  
      // Selected points
      if(v->selected) {
        pselVA.push_back(v->saveId);
        selectV.push_back(v);
      }
 
      // Border points
      if(v->margin) {
        pbrdVA.push_back(v->saveId);
        pcellVA.push_back(v->saveId);
      } else if(v->label == -1)
        pcellVA.push_back(v->saveId);
 
      // Make arrays for lines
      forall(const vertex& n, S.neighbors(v)) {
        if(S.uniqueLine(v, n)) {
          // Border lines
          if(lcnt + 1 >= lineVA.size()) {
            return;
          } 
          lineVA[lcnt++] = v->saveId;
          lineVA[lcnt++] = n->saveId;
          if(v->margin and n->margin) {
            lbrdVA.push_back(v->saveId);
            lbrdVA.push_back(n->saveId);
            lcellVA.push_back(v->saveId);
  					lcellVA.push_back(n->saveId);
          } else if(v->label == -1 and n->label == -1) {
            lcellVA.push_back(v->saveId);
            lcellVA.push_back(n->saveId);
          }
          if(v->selected or n->selected) {
            lselVA.push_back(v->saveId);
            lselVA.push_back(n->saveId);
          }
        }
      }
    }
    // Write color to GPU, the rest is done with glDrawElements
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, pcolVAid);
    glfuncs->glBufferData(GL_ARRAY_BUFFER_ARB, pcolVA.size() * sizeof(Point3GLub), &pcolVA[0], 
                                                                      GL_STATIC_DRAW_ARB); 
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);  

    REPORT_GL_ERROR("updateLines");
  }
  
  void ImgData::updateLineColor()
  {
    const vvGraph &S = mesh->graph();
    pcolVA.resize(S.size());
    if(S.empty())
      return;

    #pragma omp parallel for
    for(uint i = 0; i < S.size(); i++) {
      vertex v = S[i];
      pcolVA[i] = pColor(v);
    }

    // write to GPU
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, pcolVAid);
    glfuncs->glBufferData(GL_ARRAY_BUFFER_ARB, pcolVA.size() * sizeof(Point3GLub), &pcolVA[0], 
                                                                            GL_STATIC_DRAW_ARB); 
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);

    REPORT_GL_ERROR("updateLineColor");
  }

  // Triangle list used for selection, each gets a unique color
  void ImgData::updateSelectTris()
  {
    if(idVA.empty())
      return;

    if(idVA.size() != selVA.size()) {
      Information::out 
           << "updateSelectTris::Error: point count mismatch in sel vertex array" << endl;
      return;
    }  
  
    uint triangles = idVA.size()/3;
    #pragma omp parallel for
    for(uint i = 0; i < triangles; i++) {
      uint tcnt = i * 3;
      selVA[tcnt++] = vMapColor(i);
      selVA[tcnt++] = vMapColor(i);
      selVA[tcnt] = vMapColor(i);
    }
    // write to GPU
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, selVAid);
    glfuncs->glBufferData(GL_ARRAY_BUFFER_ARB, selVA.size() * sizeof(Point3GLub), &selVA[0], 
                                                                          GL_STATIC_DRAW_ARB); 
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);

    REPORT_GL_ERROR("updateSelectTris");
  }

  void ImgData::updateImageTex2d()
  {
    if(idVA.empty())
      return;

    if(idVA.size() != imgVA.size()) {
      Information::out 
              << "updateImageTex2d::Error: point count mismatch in sel vertex array" << endl;
      return;
    }  

    const VtxPoint2fAttr &txPos = mesh->texCoord2d();
    if(txPos.size() > 0) {
      uint triangles = idVA.size()/3;
      #pragma omp parallel for
      for(uint i = 0; i < triangles; ++i) {
        uint tcnt = i * 3;
        vertex v(idVA[tcnt]), n(idVA[tcnt + 1]), m(idVA[tcnt + 2]);
        imgVA[tcnt++] = txPos[v];
        imgVA[tcnt++] = txPos[n];
        imgVA[tcnt] = txPos[m];
      }
    }

    // write to GPU
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, imgVAid);
    glfuncs->glBufferData(GL_ARRAY_BUFFER_ARB, imgVA.size() * sizeof(Point2f), 
                                                             &imgVA[0], GL_STATIC_DRAW_ARB); 
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);

    REPORT_GL_ERROR("updateImageTex2d");
  }
   
  void ImgData::updateCellGraph()
  { 
    cellGraphPosVA.clear();
    cellGraphLineVA.clear();
    cellGraphPointVA.clear();
    cellGraphColorVA.clear();
    cellGraph &C = mesh->cells();
    cellGraphPosVA.resize(C.size());
    cellGraphPointVA.resize(C.size());
    cellGraphColorVA.resize(C.size());
    bool Mesh3D = mesh->mesh3D();
    if(C.empty())
      return;
    for(uint i = 0; i < C.size(); ++i) {
      cell c = C[i];
      c->saveId = i;
      cellGraphPosVA[i] = Point3f(c->pos - SurfOffset * c->nrml);
      cellGraphPointVA[i] = c->saveId;
      if(Mesh3D)
        cellGraphColorVA[i] = pColor(c->S.any());
      else
        cellGraphColorVA[i] = pColor(mesh->tissue().getVtx(c));
    }
    forall(const cell &c, C) {
      forall(const cell &n, C.neighbors(c)) {
        if(!C.uniqueLine(c, n))
          continue;
        cellGraphLineVA.push_back(c->saveId);
        cellGraphLineVA.push_back(n->saveId);
      }
    }   

    // Write cell positions and colors to GPU, the rest is done with glDrawElements
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, cellGraphPosVAid);
    glfuncs->glBufferData(GL_ARRAY_BUFFER_ARB, 
        cellGraphPosVA.size() * sizeof(Point3f), &cellGraphPosVA[0], GL_STATIC_DRAW_ARB); 
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);  

    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, cellGraphColorVAid);
    glfuncs->glBufferData(GL_ARRAY_BUFFER_ARB, 
        cellGraphColorVA.size() * sizeof(Point3GLub), &cellGraphColorVA[0], GL_STATIC_DRAW_ARB); 
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);

    REPORT_GL_ERROR("updateCellGraph");
  }

  void ImgData::updateAxis()
  { 
    axisPosVA.clear();
    axisColorVA.clear();

    // If there is a color array use it
    const IntVec3ColorbAttr &cellAxisColor = mesh->cellAxisColor();

    if(mesh->axisView() == "Cell Axis") {
      // Check if there is cell axis visualization data
      const IntMatrix3fAttr &cellAxisVis = mesh->cellAxisVis();
      if(cellAxisVis.empty())
        return;

      // Find the labels to draw
      IntSet labels;
      //if(mesh->hasCellGraph())
      //  forall(const cell &c, mesh->tissue().C)
      //    labels.insert(c->label)
      //else
        forall(const vertex &v, mesh->graph())
          if(v->label > 0)
            labels.insert(mesh->getLabel(v->label, mesh->parents()));

      float axisOffset = mesh->axisOffset();
  
      // Get the centers and normals
      const IntPoint3fAttr &centers
					  = mesh->useParents() ? mesh->parentCenterVis() : mesh->labelCenterVis();
      const IntPoint3fAttr &normals
            = mesh->useParents() ? mesh->parentNormalVis() : mesh->labelNormalVis();

      // Draw the 3 axes for each cell
      forall(const int label, labels) {
        const Matrix3f &axis = cellAxisVis[label];
        const Point3f &center = centers[label];
        const Point3f &normal = normals[label];
        const Vec3Colorb &colors = cellAxisColor[label];

        Point3f offset = normal * axisOffset;
  
        for(size_t i = 0; i < 3; ++i) {
          // Don't draw if axis 0
          Point3f ax = axis[i];
          if(norm(ax) == 0)
            continue;

          Point3GLub color = Point3GLub(colors[i]);
          axisColorVA.push_back(color);
          axisColorVA.push_back(color);
          axisColorVA.push_back(color);
          axisColorVA.push_back(color);
  
          // normalize the offset to make sure the angle will be the same for all axis
          Point3f off = offset * norm(ax);
          axisPosVA.push_back(center);
          axisPosVA.push_back(center - ax + off);
          axisPosVA.push_back(center);
          axisPosVA.push_back(center + ax + off);
        }
      }
    } else if(mesh->axisView() == "Triangle Axis") {
      const TriMatrix3fAttr &triangleAxisVis = mesh->triangleAxisVis();
      // Check if there is triangle visualization data
      if(triangleAxisVis.size() == 0)
        return;
      const TriVec3ColorbAttr &axisColor = mesh->triangleAxisColor();
      float axisOffset = mesh->axisOffset();

      // Loop over the triangles
      if(mesh->meshType() == "MGX3D") {
        const cellGraph &Cells = mesh->cells();
        forall(const cell &c, Cells)
          forall(const vertex &v, c->S)
            forall(const vertex &n, c->S.neighbors(v)) {
              vertex m = c->S.nextTo(v, n);
              if(!c->S.uniqueTri(v, n ,m))
                continue;
     
              Triangle t(v,n,m);
              const Matrix3f &axis = triangleAxisVis[t];
    
              Point3f center((v->pos + n->pos + m->pos)/3.0);
              Point3f normal(normalized((n->pos - v->pos) ^ (m->pos - v->pos)));
    
              const Vec3Colorb &colors = axisColor[t];
    
              Point3f offset = normal * axisOffset;
      
              for(size_t i = 0; i < 3; ++i) {
                // Don't draw if axis 0 or aplha 0
                Point3f ax = axis[i];
                if(norm(ax) == 0) // or colors[i][3] == 0)
                  continue;
    
                Point3GLub color = Point3GLub(colors[i]);
                axisColorVA.push_back(color);
                axisColorVA.push_back(color);
                axisColorVA.push_back(color);
                axisColorVA.push_back(color);
        
                // normalize the offset to make sure the angle will be the same for all axis
                Point3f off = offset * norm(ax);
                axisPosVA.push_back(center);
                axisPosVA.push_back(center - ax + off);
                axisPosVA.push_back(center);
                axisPosVA.push_back(center + ax + off);
              }
            }
      } else {
        const vvGraph &S = mesh->graph();
        forall(const vertex &v, S)
          forall(const vertex &n, S.neighbors(v)) {
            vertex m = S.nextTo(v, n);
            if(!S.uniqueTri(v, n ,m))
              continue;
   
            Triangle t(v,n,m);
            const Matrix3f &axis = triangleAxisVis[t];
  
            Point3f center((v->pos + n->pos + m->pos)/3.0);
            Point3f normal(normalized((n->pos - v->pos) ^ (m->pos - v->pos)));
  
            const Vec3Colorb &colors = axisColor[t];
  
            Point3f offset = normal * axisOffset;
    
            for(size_t i = 0; i < 3; ++i) {
              // Don't draw if axis 0 or aplha 0
              Point3f ax = axis[i];
              if(norm(ax) == 0) // or colors[i][3] == 0)
                continue;
  
              Point3GLub color = Point3GLub(colors[i]);
              axisColorVA.push_back(color);
              axisColorVA.push_back(color);
              axisColorVA.push_back(color);
              axisColorVA.push_back(color);
      
              // normalize the offset to make sure the angle will be the same for all axis
              Point3f off = offset * norm(ax);
              axisPosVA.push_back(center);
              axisPosVA.push_back(center - ax + off);
              axisPosVA.push_back(center);
              axisPosVA.push_back(center + ax + off);
            }
          }
      }
    } else if(mesh->axisView() == "Vertex Axis") {
      // Check if there is cell axis visualization data
      const VtxMatrix3fAttr &vertexAxisVis = mesh->vertexAxisVis();
      if(vertexAxisVis.empty())
        return;

      float axisOffset = mesh->axisOffset();
      const VtxVec3ColorbAttr &axisColor = mesh->vertexAxisColor();
			
      // Draw the 3 axes for each vertex
			const vvGraph &S = mesh->graph();
      forall(const vertex &v, S) {
        const Matrix3f &axis = vertexAxisVis[v];
        const Point3f center(v->pos);
        const Point3f normal(v->nrml);
        const Vec3Colorb &colors = axisColor[v];

        Point3f offset = normal * axisOffset;
  
        for(size_t i = 0; i < 3; ++i) {
          // Don't draw if axis 0
          Point3f ax = axis[i];
          if(norm(ax) == 0)
            continue;

          Point3GLub color = Point3GLub(colors[i]);
          axisColorVA.push_back(color);
          axisColorVA.push_back(color);
          axisColorVA.push_back(color);
          axisColorVA.push_back(color);
  
          // normalize the offset to make sure the angle will be the same for all axis
          Point3f off = offset * norm(ax);
          axisPosVA.push_back(center);
          axisPosVA.push_back(center - ax + off);
          axisPosVA.push_back(center);
          axisPosVA.push_back(center + ax + off);
        }
      }
    }

    // Write axis positions and colors to GPU
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, axisPosVAid);
    glfuncs->glBufferData(GL_ARRAY_BUFFER_ARB, axisPosVA.size() * sizeof(Point3f), &axisPosVA[0], GL_STATIC_DRAW_ARB); 
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);  

    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, axisColorVAid);
    glfuncs->glBufferData(GL_ARRAY_BUFFER_ARB, axisColorVA.size() * sizeof(Point3GLub), &axisColorVA[0], GL_STATIC_DRAW_ARB); 
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);

    REPORT_GL_ERROR("updateAxis");
  }


  // Setup framebuffer for selection
  bool ImgData::startSelectFbo()
  {
    // Size of renderbuffers
    static GLint w = 0, h = 0;

    // Get viewport size
    GLint dims[4] = {0};
    glGetIntegerv(GL_VIEWPORT, dims);
    GLint width = dims[2];
    GLint height = dims[3];

    // If there is no frame buffer create one
    if(!selFboId)
      glfuncs->glGenFramebuffers(1, &selFboId);
    glfuncs->glBindFramebuffer(GL_DRAW_FRAMEBUFFER, selFboId);

    // (re-)create if not created or changed in size
    bool recreate = w != width or h != height or !selFboColor or !selFboDepth;
    if(recreate) {
      if(selFboColor)
        glfuncs->glDeleteRenderbuffers(1, &selFboColor);

      glfuncs->glGenRenderbuffers(1, &selFboColor);
      glfuncs->glBindRenderbuffer(GL_RENDERBUFFER_EXT, selFboColor);
      glfuncs->glRenderbufferStorage(GL_RENDERBUFFER_EXT, GL_RGB8, width, height);

      if(selFboDepth)
        glfuncs->glDeleteRenderbuffers(1, &selFboDepth);
      glfuncs->glGenRenderbuffers(1, &selFboDepth);
      glfuncs->glBindRenderbuffer(GL_RENDERBUFFER, selFboDepth);
      glfuncs->glRenderbufferStorage(GL_RENDERBUFFER, GL_DEPTH_COMPONENT, width, height);
      glfuncs->glFramebufferRenderbuffer(GL_DRAW_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_RENDERBUFFER, selFboColor);
      glfuncs->glFramebufferRenderbuffer(GL_DRAW_FRAMEBUFFER, GL_DEPTH_ATTACHMENT, GL_RENDERBUFFER, selFboDepth);
      // Check frame buffer created properly
      if(glfuncs->glCheckFramebufferStatus(GL_DRAW_FRAMEBUFFER) != GL_FRAMEBUFFER_COMPLETE) {
        Information::out << "Warning: Unable to create framebuffer, selection will not work." << endl;
        return false;
      }
      glViewport(0, 0, width, height);
      w = width;
      h = height;
    }

    REPORT_GL_ERROR("startSelectFbo");

    return true;
  }

  bool ImgData::stopSelectFbo()
  {
    glfuncs->glBindFramebuffer(GL_DRAW_FRAMEBUFFER, 0);
    glfuncs->glBindFramebuffer(GL_READ_FRAMEBUFFER, 0);

    REPORT_GL_ERROR("stopSelectFbo");

    return true;
  }

  bool ImgData::readSelectFbo()
  {
    if(!selFboId)
       startSelectFbo();

    glfuncs->glBindFramebuffer(GL_READ_FRAMEBUFFER, selFboId);
    glReadBuffer( GL_COLOR_ATTACHMENT0 );

    REPORT_GL_ERROR("readSelectFboUnbind");

    return true;
  }

  // Fill the vertex array
  void ImgData::fillVBOs()
  {  
    // If nothing to draw clear arrays and exit
    if(mesh->empty()) {
      idVA.clear();
      triCell.clear();

      posVA.clear();
      nrmlVA.clear();
      selVA.clear();
      texVA.clear();
      imgVA.clear();
      triangleColorVA.clear();
      cellColorVA.clear();

      pntsVA.clear();
      pcolVA.clear();
      pbrdVA.clear();
      pcellVA.clear();
      pselVA.clear();
  
      lineVA.clear();
      lbrdVA.clear();
      lcellVA.clear();
      lselVA.clear();

      cellGraphPosVA.clear();
      cellGraphColorVA.clear();
      cellGraphPointVA.clear();
      cellGraphLineVA.clear();

      axisPosVA.clear();
      axisColorVA.clear();

      return;
    }

#ifdef DEBUG_CC_TIMING
    QTime vboTimer;
    vboTimer.start(); 
#endif

    try {

      // Set up Vertex Buffer Objects to put data onto graphics card
      if(!posVAid)
        glfuncs->glGenBuffers(1, &posVAid);
      if(!nrmlVAid)
        glfuncs->glGenBuffers(1, &nrmlVAid);
      if(!selVAid)
        glfuncs->glGenBuffers(1, &selVAid);
      if(!texVAid)
        glfuncs->glGenBuffers(1, &texVAid);
      if(!imgVAid)
        glfuncs->glGenBuffers(1, &imgVAid);
      if(!triangleColorVAid)
        glfuncs->glGenBuffers(1, &triangleColorVAid);
      if(!cellColorVAid)
        glfuncs->glGenBuffers(1, &cellColorVAid);

      if(!pntsVAid)
        glfuncs->glGenBuffers(1, &pntsVAid);
      if(!pcolVAid)
        glfuncs->glGenBuffers(1, &pcolVAid);
      if(!pselVAid)
        glfuncs->glGenBuffers(1, &pselVAid);  

      if(!cellGraphPosVAid)
        glfuncs->glGenBuffers(1, &cellGraphPosVAid);
      if(!cellGraphColorVAid)
        glfuncs->glGenBuffers(1, &cellGraphColorVAid);

      if(!axisPosVAid)
        glfuncs->glGenBuffers(1, &axisPosVAid);
      if(!axisColorVAid)
        glfuncs->glGenBuffers(1, &axisColorVAid);

      // Get the graph
      vvGraph &S = mesh->graph();

      // Also counts triangles and fills in idVA[]
      uint VCount = 0, LCount = 0, TCount = 0;
      updatePoints(VCount, LCount, TCount);

      // Set the normals and mark the margins, no point for 3D
      bool mesh3D = mesh->mesh3D();
      if(!mesh3D) {
        setNormals(S);
        markMargins(S);
      }

      // Check the graph
      mesh->chkGraph(S);

      // Resize vertex arrays
      pcolVA.resize(VCount);
      lineVA.resize(LCount * 2);
      posVA.resize(TCount * 3);
      nrmlVA.resize(TCount * 3);
      selVA.resize(TCount * 3);
      texVA.resize(TCount * 3);
      imgVA.resize(TCount * 3);

      // Now fill the VBOs
      updateTriPos();
      updateTriColor();
      updateLines();
      updateSelectTris();
      updateImageTex2d();
      updateCellGraph();
      updateAxis();
      //mesh->complexHandler().update();

    } catch(QString msg) {
      Information::setStatus(QString("fillVBOs::Error: %1").arg(msg));
    } catch(...) {
      Information::setStatus("fillVBOs::Error unknown exception");
    }

#ifdef DEBUG_CC_TIMING
    float vboSec = vboTimer.elapsed()/1000.0;
    Information::out << "VV VBO create time:" << vboSec << " seconds" << endl;
#endif
  }

  // Find selected triangle in triangle select mode (!MeshSelect)
  void ImgData::findSelectTriangle(uint x, uint y, std::vector<uint> &vList, int &label, bool useParentLabel)
  {
    // Draw unique colors in the select frame buffer
    drawSurf(true);

    // Bind select framebuffer for reading
    readSelectFbo();
  
    // Find the color of the pixels drawn in BACK buffer
    Point3GLub pix;
    GLint viewport[4];
    glGetIntegerv(GL_VIEWPORT, viewport);
    glReadPixels(x, viewport[3] - y, 1, 1, GL_RGB, GL_UNSIGNED_BYTE, &pix);

    // Unbind select framebuffer
    stopSelectFbo();

    int u = vMapColor(pix);

    // Return if nothing selected
    if(u < 0)
      return;
  
    u *= 3;
    if(u + 3 > int(idVA.size())) {
      Information::out << "findSelectTriangle:Error Index too large, tri idx:" 
                       << u << " idVA size:" << idVA.size() << endl;
      return;
    }
    // Return triangle and label for selected color
    for(uint i = 0; i < 3; i++)
      vList.push_back(u + i);
    if(useParentLabel)
      label = mesh->getLabel(idVA[u], idVA[u + 1], idVA[u + 2], mesh->parents());
    else
      label = getLabel(idVA[u], idVA[u + 1], idVA[u + 2]);
  }
  
  // Trim bounding box with a clip plane in point normal form
  void ImgData::bBoxClip(BoundingBox3f& bBox, Point3f p, Point3f n)
  {
    int x, y, z;
    Point3f u1, u2, pt;
    float dir, s;
  
    // check x direction, find which line to test
    dir = n * Point3f(1.0f, 0.0f, 0.0f);
    y = (n.y() > 0 ? 1 : 0);
    z = (n.z() > 0 ? 1 : 0);
    u1 = Point3f(bBox[0].x(), bBox[y].y(), bBox[z].z());
    u2 = Point3f(bBox[1].x(), bBox[y].y(), bBox[z].z());
    if(planeLineIntersect(p, n, u1, u2, s, pt)) {
      if(dir > 0) {     // outside points in positive x dir
        if(pt.x() > bBox[0].x())
          bBox[0].x() = pt.x();
      } else {
        if(pt.x() < bBox[1].x())
          bBox[1].x() = pt.x();
      }
    }
    // check y direction, find which line to test
    dir = n * Point3f(0.0f, 1.0f, 0.0f);
    x = (n.x() > 0 ? 1 : 0);
    z = (n.z() > 0 ? 1 : 0);
    u1 = Point3f(bBox[x].x(), bBox[0].y(), bBox[z].z());
    u2 = Point3f(bBox[x].x(), bBox[1].y(), bBox[z].z());
    if(planeLineIntersect(p, n, u1, u2, s, pt)) {
      if(dir > 0) {
        if(pt.y() > bBox[0].y())
          bBox[0].y() = pt.y();
      } else {
        if(pt.y() < bBox[1].y())
          bBox[1].y() = pt.y();
      }
    }
    // check z direction, find which line to test
    dir = n * Point3f(0.0f, 0.0f, 1.0f);
    x = (n.x() > 0 ? 1 : 0);
    y = (n.y() > 0 ? 1 : 0);
    u1 = Point3f(bBox[x].x(), bBox[y].y(), bBox[0].z());
    u2 = Point3f(bBox[x].x(), bBox[y].y(), bBox[1].z());
    if(planeLineIntersect(p, n, u1, u2, s, pt)) {
      if(dir > 0) {
        if(pt.z() > bBox[0].z())
          bBox[0].z() = pt.z();
      } else {
        if(pt.z() < bBox[1].z())
          bBox[1].z() = pt.z();
      }
    }
  }
  
  //  Calculate bounding box from clipping planes (in point normal form)
  void ImgData::bBoxFromClip()
  {
    // clear, set to entire texture
    Point3f origin = stack->origin();
    const Point3f size(stack->size());
    bBox[0] = origin;
    bBox[1] = origin + size;
  
    // Test all 6 planes (twice)
    for(int i = 0; i < 2; i++)
      for(int pln = 0; pln < 6; pln++) {
        if(!clipDo[pln / 2])
          continue;
  
        Point3f n(pn[pln].x(), pn[pln].y(), pn[pln].z());
        float d = pn[pln][3];
        bBoxClip(bBox, n * -d, n);
      }
  }
  
  // Setup data for clipping test
  void ImgData::getClipTestData(ClipRegion& clip1, ClipRegion& clip2, ClipRegion& clip3)
  {
    Clip* c1 = clip1.clip;
    Clip* c2 = clip2.clip;
    Clip* c3 = clip3.clip;
    // Indicator if in use
    clipDo[0] = (c1->enabled() or c1->grid() ? 1 : 0);
    clipDo[1] = (c2->enabled() or c2->grid() ? 1 : 0);
    clipDo[2] = (c3->enabled() or c3->grid() ? 1 : 0);
  
    // Frame matrix
    frm = Matrix4d(getFrame().worldMatrix());
  
    // Clip plane matrices
    clm[0] = Matrix4d(c1->frame().inverse().worldMatrix());
    clm[1] = Matrix4d(c2->frame().inverse().worldMatrix());
    clm[2] = Matrix4d(c3->frame().inverse().worldMatrix());
  
    // Clipping planes, point normal form
    pn[0] = Point4f(frm * clm[0] * 
                 Point4d(c1->normal().x(), c1->normal().y(), c1->normal().z(), c1->width()));
    pn[1] = Point4f(frm * clm[0] * 
                 Point4d(-c1->normal().x(), -c1->normal().y(), -c1->normal().z(), c1->width()));
    pn[2] = Point4f(frm * clm[1] * 
                 Point4d(c2->normal().x(), c2->normal().y(), c2->normal().z(), c2->width()));
    pn[3] = Point4f(frm * clm[1] * 
                 Point4d(-c2->normal().x(), -c2->normal().y(), -c2->normal().z(), c2->width()));
    pn[4] = Point4f(frm * clm[2] * 
                 Point4d(c3->normal().x(), c3->normal().y(), c3->normal().z(), c3->width()));
    pn[5] = Point4f(frm * clm[2] * 
                 Point4d(-c3->normal().x(), -c3->normal().y(), -c3->normal().z(), c3->width()));
  
    // Put in host vector
    Hpn.resize(6);
    for(int i = 0; i < 6; i++)
      Hpn[i] = pn[i];
  }
  
  void ImgData::voxelEditStart(ClipRegion& clip1, ClipRegion& clip2, ClipRegion& clip3)
  {
    if(stack->work()->isVisible()) {
      // Setup matrices and planes for clip test
      getClipTestData(clip1, clip2, clip3);
  
      // If 3D seeding, increment label number
      if(FillWorkData and stack->work()->isVisible() and SeedStack and stack->work()->labels()) {
        int lab = stack->nextLabel();
        Information::setStatus(QString("Adding seed %1").arg(lab));
      }
  
      // Setup texture update area
      bBoxTex[0] = Point3i(stack->size().x(), stack->size().y(), stack->size().z());
      bBoxTex[1] = Point3i(0, 0, 0);
    }
  }
 
  // Edit the 3D voxel data 
  void ImgData::voxelEditStop()
  {
    // Reload the texture to see the edit result
    if(stack->work()->isVisible() and pixelsChanged)
      updateTex(workDataTexId, bBoxTex);
    pixelsChanged = false;
  }
  
  void ImgData::voxelEdit(float pixelRadius, const Point3f& p, const Point3f& px, 
                       const Point3f& py, const Point3f& pz, bool doCut, int currentLabel)
  {
    Point3f Origin = stack->origin();
    if(stack->work()->isVisible()) {
      // Clear bounding box, set to entire texture
      bBox[0] = Origin;
      bBox[1] = Origin + multiply(Point3f(stack->size()), stack->step());
      // Trim bounding box to clip planes
      for(int i = 0; i < 2; i++) {
        // Test all 6 planes (twice)
        for(int pln = 0; pln < 6; pln++) {
          if(!clipDo[pln / 2])
            continue;
  
          Point3f n(pn[pln].x(), pn[pln].y(), pn[pln].z());
          float d = pn[pln][3];
          bBoxClip(bBox, n * -d, n);
        }
        bBoxClip(bBox, p - pixelRadius * px, px);
        bBoxClip(bBox, p + pixelRadius * px, -px);
        bBoxClip(bBox, p - pixelRadius * py, py);
        bBoxClip(bBox, p + pixelRadius * py, -py);
        if(doCut) {
          bBoxClip(bBox, p - pixelRadius * pz, pz);
          bBoxClip(bBox, p + pixelRadius * pz, -pz);
        }
      }
  
      // Get bounding box in image coordinates and pad a pixel
      Point3i imgSize(stack->size().x(), stack->size().y(), stack->size().z());
      Point3i base = worldToImagei(bBox[0]) - Point3i(1, 1, 1);
      Point3i end = worldToImagei(bBox[1]) + Point3i(1, 1, 1);
  
      // Check bounds and update bound box for texture update
      for(int i = 0; i < 3; i++) {
        base[i] = std::max(base[i], 0);
        end[i] = std::min(end[i], imgSize[i]);
        bBoxTex[0][i] = std::min(bBoxTex[0][i], base[i]);
        bBoxTex[1][i] = std::max(bBoxTex[1][i], end[i]);
      }
      Point3i size = end - base;
  
      // Call cuda
      if(size.x() > 0 and size.y() > 0 and size.z() > 0) {
        ushort pixval = 0x0;
        if(FillWorkData) {
          if(stack->work()->labels()) {
            if(SeedStack)
              pixval = stack->viewLabel();
            else
              pixval = currentLabel;
          } else
            pixval = 0xFFFF;
        }
        mgx::voxelEditGPU(base, size, imgSize, stack->step(), Origin, p, pz, pixelRadius, 
                                     pixval, clipDo, Hpn, stack->work()->data());
        pixelsChanged = true;
      }
      if(pixelsChanged
         and size_t(bBoxTex[1].x() - bBoxTex[0].x()) * (bBoxTex[1].y() - bBoxTex[0].y())
         * (bBoxTex[1].z() - bBoxTex[0].z()) < size_t(VoxelEditMaxPix)) {
        updateTex(workDataTexId, bBoxTex);
        pixelsChanged = false;
        bBoxTex[0] = Point3i(stack->size().x(), stack->size().y(), stack->size().z());
        bBoxTex[1] = Point3i(0, 0, 0);
      }
    }
  }
  
	// Find the seed point for 3D voxel editing
	bool ImgData::findSeedPoint(uint x, uint y, CutSurf& cutSurf, Point3f& p)
  {
    if(!SeedStack)
      return (false); 
  
    vvGraph& S = mesh->graph();
    p = Point3f(0, 0, 0);
    // Use the cutting surface if visible to locate the seed in Z
    if(cutSurf.cut->isVisible()) {
      // Bind select frame buffer for reading, do we need to clear?
      readSelectFbo();
      cutSurf.drawCutSurf(*this, true);

      // Pixel color is x,y,z
      GLint viewport[4];
      glGetIntegerv(GL_VIEWPORT, viewport);
      glReadPixels(x, viewport[3] - y, 1, 1, GL_RGB, GL_FLOAT, &p);

      // Unbind select frame buffer
      stopSelectFbo();
      p -= Point3f(.5f, .5f, .5f);
    } else if(mesh->showSurface() and !S.empty()) {
      // Otherwise use the mesh surface if enabled
      std::vector<uint> vList;
      int label;
      findSelectTriangle(x, y, vList, label);
      if(vList.empty())
        return (false);
      forall(uint u, vList) {
        if(u < idVA.size()) {
          vertex v = idVA[u];
          p += Point3f(v->pos);
        }
      }
      p /= vList.size();
    }
    return (true);
  }
  
  void ImgData::updateSelection()
  {
    selectV.clear();
    pselVA.clear();
    lselVA.clear();

    vvGraph &S = mesh->graph();

    if(S.empty())
      return;

    forall(const vertex &v, S) {
      if(v->selected) {
        selectV.push_back(v);
        pselVA.push_back(v->saveId);
        forall(const vertex& n, S.neighbors(v)) {
          lselVA.push_back(v->saveId);
          lselVA.push_back(n->saveId);
        }
      }
    }
    updateLineColor();
    updateCellGraph();
  }
  
  // Clear mesh selection
  void ImgData::clearMeshSelect()
  {
    if(mesh->clearSelection())
      updateSelection();
  }
  
  // Remove points from selection
  void ImgData::removeSelect(const VtxVec& vList)
  {
    bool changed = false;
    #pragma omp parallel for
    for(uint i = 0; i < vList.size(); i++) {
      vertex v = vList[i];
      if(v->selected) {
        v->selected = false;
        changed = true;
      }
    }
    if(changed)
      updateSelection();
  }
  
  // Add points to selection
  void ImgData::addSelect(const VtxVec& vList)
  {
    bool changed = false;
    #pragma omp parallel for
    for(uint i = 0; i < vList.size(); i++) {
      vertex v = vList[i];
      if(!v->selected) {
        v->selected = true;
        changed = true;
      }
    }
    if(changed)
      updateSelection();

    // In special cases print some useful info
    if(selectV.size() == 1)
      Information::setStatus(QString("Vertex label: %1").arg(selectV[0]->label));
    else if(selectV.size() == 2) {
      vertex v1 = selectV[0];
      vertex v2 = selectV[1];
      Information::setStatus(QString("Distance between vertices: %1um").arg((v1->pos - v2->pos).norm())); 
    } else
      Information::setStatus(QString("%1 vertex(es) selected").arg(selectV.size()));
  }

  // Update a list of vertex labels
  void ImgData::updLabel(int label, const std::vector<uint> &vList)
  {
    // Update vertex arrays
    bool changed = false;
    vvGraph M;
    IntFloatAttr &labelHeat = mesh->labelHeat();
    bool isCellMesh = (mesh->meshType() == "MGXC");
    forall(uint u, vList) {
      vertex v(idVA[u]);
      if(!isCellMesh or v->type == 'c') {
        v->label = label;
        M.insert(v);
        changed = true;
      }
      texVA[u] = texCoordLabels(label, v, labelHeat);
    }
    if(!changed)
      return;
  
    markMargins(M, mesh->graph());
  
    // Write texture data to video memory
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, texVAid);

    // Try using mapped memory to just write what we need to
    Point3f* tp = (Point3f*)(glfuncs->glMapBuffer(GL_ARRAY_BUFFER_ARB, GL_WRITE_ONLY_ARB));
    if(!tp) {
      REPORT_GL_ERROR("glMapBufferARB(texVAid) - using glBufferData");
      glfuncs->glBufferData(GL_ARRAY_BUFFER_ARB, texVA.size() * sizeof(Point3f), &texVA[0], 
                                                                 GL_STATIC_DRAW_ARB);
    } else {
      forall(uint u, vList)
        tp[u] = texVA[u];
      glfuncs->glUnmapBuffer(GL_ARRAY_BUFFER_ARB);
    }
  
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);
  }

  // Update a list of vertex labels
  void ImgData::updSignal(int signal, const std::vector<uint> &vList)
  {
    // Update vertex arrays
    IntFloatAttr &labelHeat = mesh->labelHeat();
    forall(uint u, vList) {
      vertex v(idVA[u]);
      v->signal = signal;
      changed = true;
      texVA[u] = texCoordLabels(v->label, v, labelHeat);
    }

    // Write texture data to video memory
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, texVAid);

    // Try using mapped memory to just write what we need to
    Point3f* tp = (Point3f*)(glfuncs->glMapBuffer(GL_ARRAY_BUFFER_ARB, GL_WRITE_ONLY_ARB));
    if(!tp) {
      REPORT_GL_ERROR("glMapBufferARB(texVAid) - using glBufferData");
      glfuncs->glBufferData(GL_ARRAY_BUFFER_ARB, texVA.size() * sizeof(Point3f), &texVA[0], 
                                                                 GL_STATIC_DRAW_ARB);
    } else {
      forall(uint u, vList)
        tp[u] = texVA[u];
      glfuncs->glUnmapBuffer(GL_ARRAY_BUFFER_ARB);
    }
  
    glfuncs->glBindBuffer(GL_ARRAY_BUFFER_ARB, 0); 
  }

  void ImgData::updateHistogram(std::vector<double>& hist, const HVecUS& data, 
                          std::pair<double, double>& minMaxValues, int max_data, int size)
  {
    hist.clear();
    if(data.empty())
      return;
    hist.resize(size, 0);
    uint min = 65536;
    uint max = 0;
    #pragma omp parallel for reduction(min:min) reduction(max:max)
    for(uint i = 0; i < data.size(); i++) {
      ushort d = data[i];
      int val = (d * size) / max_data;
      hist[val]++;
      if(min > d)
        min = d;
      if(max < d)
        max = d;
    }
    minMaxValues.first = min / double(max_data);
    minMaxValues.second = max / double(max_data);
  }
  
  void ImgData::updateSurfHistogram()
  {
    if(SurfHist.empty()) {
      // TODO
    }
  }
  
  void ImgData::updateHeatHistogram()
  {
    if(HeatHist.empty()) {
      // TODO
    }
  }
  void ImgData::updateWorkHistogram()
  {
    if(WorkHist.empty())
      updateHistogram(WorkHist, stack->work()->data(), workBounds);
  }
  
  void ImgData::updateMainHistogram()
  {
    if(MainHist.empty())
      updateHistogram(MainHist, stack->main()->data(), mainBounds);
  }
  
  void ImgData::editMainTransferFunction()
  {
    updateMainHistogram();
    mainTransferDlg->changeHistogram(MainHist);
    mainTransferDlg->changeBounds(mainBounds);
    mainTransferDlg->setTransferFunction(stack->main()->transferFct());
    if(mainTransferDlg->exec() == QDialog::Accepted) {
      emit changedInterface();
      stack->main()->setTransferFct(mainTransferDlg->transferFunction());
      updateMainColorMap();
    }
  }
  
  void ImgData::editWorkTransferFunction()
  {
    updateWorkHistogram();
    workTransferDlg->changeHistogram(WorkHist);
    workTransferDlg->changeBounds(workBounds);
    workTransferDlg->setTransferFunction(stack->work()->transferFct());
    if(workTransferDlg->exec() == QDialog::Accepted) {
      emit changedInterface();
      stack->work()->setTransferFct(workTransferDlg->transferFunction());
      updateWorkColorMap();
    }
  }
   
  void ImgData::editSurfTransferFunction()
  {
    updateSurfHistogram();
    surfTransferDlg->changeHistogram(SurfHist);
    surfTransferDlg->setTransferFunction(mesh->surfFct());
    if(surfTransferDlg->exec() == QDialog::Accepted) {
      emit changedInterface();
      mesh->setSurfFct(surfTransferDlg->transferFunction());
      updateSurfColorMap();
    }
  }
  
  void ImgData::editHeatTransferFunction()
  {
    updateHeatHistogram();
    heatTransferDlg->changeHistogram(HeatHist);
    heatTransferDlg->setTransferFunction(mesh->heatFct());
    if(heatTransferDlg->exec() == QDialog::Accepted) {
      emit changedInterface();
      mesh->setHeatFct(heatTransferDlg->transferFunction());
      updateHeatColorMap();
    }
  } 

  void ImgData::updateColorMap(bool& newColorMap, std::vector<TransferFunction::Colorf>& ColorMap,
                               const TransferFunction& transferFct)
  {
    static const int n = 4096;
    ColorMap.resize(n);
    double dc = 1.0 / (n - 1);
    for(int i = 0; i < n; ++i) {
      ColorMap[i] = transferFct.rgba(i * dc);
    }
    newColorMap = true;
  }

  bool ImgData::setupTexFromColorMap(GLuint &texId, const ColorbVec &colorVec)
  {
    if(!texId)
      glGenTextures(1, &texId);
    glBindTexture(GL_TEXTURE_1D, texId);
    glfuncs->glTexImage1D(GL_TEXTURE_1D, 0, GL_RGBA8, colorVec.size(), 0, GL_RGBA, GL_UNSIGNED_BYTE, &colorVec[0][0]);
    glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
    glBindTexture(GL_TEXTURE_1D, 0);

    return true;
  };
  
  void ImgData::setupColorMap(bool& newColorMap, GLuint &texId, const std::vector<TransferFunction::Colorf>& ColorMap,
                              Shader::ActiveTextures activeTex)
  {
    if(newColorMap or !texId) {
      newColorMap = false;
      if(!texId)
        glGenTextures(1, &texId);
      glBindTexture(GL_TEXTURE_1D, texId);
      glfuncs->glTexImage1D(GL_TEXTURE_1D, 0, GL_RGBA16, ColorMap.size(), 0, GL_RGBA, GL_FLOAT, &ColorMap[0]);
      glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
      glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
      glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
      glBindTexture(GL_TEXTURE_1D, 0);
    }
    Shader::activeTexture(activeTex);
    glBindTexture(GL_TEXTURE_1D, texId);
    Shader::activeTexture(Shader::AT_NONE);
  }
  
  void ImgData::updateWorkColorMap(const TransferFunction& fct)
  {
    stack->work()->setTransferFct(fct);
    updateWorkColorMap();
  }
  
  void ImgData::updateMainColorMap(const TransferFunction& fct)
  {
    stack->main()->setTransferFct(fct);
    updateMainColorMap();
  }
  
  void ImgData::updateSurfColorMap(const TransferFunction& fct)
  {
    mesh->setSurfFct(fct);
    updateSurfColorMap();
  }
  
  void ImgData::updateHeatColorMap(const TransferFunction& fct)
  {
    mesh->setHeatFct(fct);
    updateHeatColorMap();
  }
  
  void ImgData::setColorMap(const QString& pth, bool work)
  {
    QFile f(pth);
    if(!f.open(QIODevice::ReadOnly)) {
      err << "Error, cannot open file '" << pth << "' for reading" << endl;
      return;
    }
    QTextStream ts(&f);
    QString content = ts.readAll();
    TransferFunction fct = TransferFunction::load(content);
    if(fct.empty()) {
      err << "Error, file '" << pth << "' doesn't contain a valid transfer function" << endl;
      return;
    }
    if(work)
      workTransferDlg->setTransferFunction(fct);
    else
      mainTransferDlg->setTransferFunction(fct);
  }
  
  void ImgData::setWork16Bit(bool val)
  {
    if(val != Work16Bit) {
      Work16Bit = val;
      reloadTex(workDataTexId);
    }
  }
  
  void ImgData::setMain16Bit(bool val)
  {
    if(val != Main16Bit) {
      Main16Bit = val;
      reloadTex(mainDataTexId);
    }
  }
  
  // Return texture info for a vertex, used when triangle context not known (messes up wall heat map)
  inline Point3f ImgData::texCoordLabels(int label, const vertex &v, const IntFloatAttr &labelHeat)
  {
    if(label <= 0)
      return Point3f(0.0f, v->signal, 0.0f);
  
    float retlab = float(label);
    // see if vertex has heat
    IntFloatAttr::const_iterator found = labelHeat.find(label);
    if(found == labelHeat.end())
      return Point3f(retlab, v->signal, 0.0f);
  
    // Vertex has a cell based heatmap
    float heat = found->second;
    // Use a negative label to represent a valid heat value
    return Point3f(-retlab, v->signal, heat);
  }
  
  // Return texture info for a vertex with triangle context. This allows to color walls.
  // Be careful adding more maps here, it is called using omp in updateTriColor
  inline Point3f ImgData::texCoordWalls(int label, const vertex &v, const vertex &a, 
                    const vertex &b, const vertex &c, const IntIntFloatAttr &wallHeat)
  {
    IntIntPair wall;
    float retlab = float(label);
    bool isBord = mesh->isBordTriangle(label, a, b, c, wall);
    // See if triangle has wall heat
    IntIntFloatAttr::const_iterator found = wallHeat.end();

    if(isBord)
      found = wallHeat.find(wall);
    if(found == wallHeat.end())
      return Point3f(retlab, v->signal, 0.0f);
  
    // Vertex has a wall heat, negative label represents valid heat
    return Point3f(-retlab, v->signal, found->second);
  }
  
  // Slots for stack controls
  // Main stack
  void ImgData::MainShowSlot(bool val)
  {
    if(val != stack->main()->isVisible()) {
      if(val)
        stack->main()->show();
      else
        stack->main()->hide();
      emit changedInterface();
      emit viewerUpdate();
    }
  }
  
  void ImgData::MainBrightSlot(int val)
  {
    float v = float(val) / 10000.0;
    if(stack->main()->brightness() != v) {
      stack->main()->setBrightness(v);
      emit changedInterface();
      emit viewerUpdate();
    }
  }
  
  void ImgData::MainOpacitySlot(int val)
  {
    float v = float(val) / 10000.0;
    if(stack->main()->opacity() != v) {
      stack->main()->setOpacity(v);
      emit changedInterface();
      emit viewerUpdate();
    }
  }
  
  void ImgData::MainLabelsSlot(bool val)
  {
    if(val != stack->main()->labels()) {
      emit changedInterface();
      stack->main()->setLabels(val);
      reloadMainTex();
      emit viewerUpdate();
    }
  }
  
  void ImgData::Main16BitSlot(bool val)
  {
    if(val != Main16Bit) {
      emit changedInterface();
      Main16Bit = val;
      reloadMainTex();
      emit viewerUpdate();
    }
  }
  
  // Work stack
  void ImgData::WorkShowSlot(bool val)
  {
    if(val != stack->work()->isVisible()) {
      emit changedInterface();
      if(val)
        stack->work()->show();
      else
        stack->work()->hide();
      emit viewerUpdate();
    }
  }
  
  void ImgData::WorkBrightSlot(int val)
  {
    float v = float(val) / 10000.0;
    if(stack->work()->brightness() != v) {
      emit changedInterface();
      stack->work()->setBrightness(v);
      emit viewerUpdate();
    }
  }
  
  void ImgData::WorkOpacitySlot(int val)
  {
    float v = float(val) / 10000.0;
    if(stack->work()->opacity() != v) {
      emit changedInterface();
      stack->work()->setOpacity(v);
      emit viewerUpdate();
    }
  }
  
  void ImgData::WorkLabelsSlot(bool val)
  {
    if(val != stack->work()->labels()) {
      emit changedInterface();
      stack->work()->setLabels(val);
      reloadWorkTex();
      emit viewerUpdate();
    }
  }
  
  void ImgData::Work16BitSlot(bool val)
  {
    if(val != Work16Bit) {
      emit changedInterface();
      Work16Bit = val;
      reloadWorkTex();
      emit viewerUpdate();
    }
  }

  void ImgData::MainColorMapSlot()
  {
    editMainTransferFunction();
    emit viewerUpdate();
  }

  void ImgData::WorkColorMapSlot()
  {
    editWorkTransferFunction();
    emit viewerUpdate();
  }

  // Surface
  void ImgData::SurfShowSlot(bool val)
  {
    if(val != mesh->showSurface()) {
      emit changedInterface();
      mesh->setShowSurface(val);
      emit viewerUpdate();
    }
  }
  
  void ImgData::SurfBrightSlot(int val)
  {
    float v = float(val) / 10000.0;
    if(mesh->brightness() != v) {
      emit changedInterface();
      mesh->setBrightness(v);
      emit viewerUpdate();
    }
  }
  
  void ImgData::SurfOpacitySlot(int val)
  {
    float v = float(val) / 10000.0;
    if(mesh->opacity() != v) {
      emit changedInterface();
      mesh->setOpacity(v);
      emit viewerUpdate();
    }
  }
  
  void ImgData::SurfBlendSlot(bool val)
  {
    if(val != mesh->blending()) {
      emit changedInterface();
      mesh->setBlending(val);
      emit viewerUpdate();
    }
  }
  
  void ImgData::SurfVertexSlot(bool val)
  {
    if(val and mesh->surfView() != Mesh::SURF_VERTEX) {
      emit changedInterface();
      mesh->setShowVertex();
      emit viewerUpdate();
    }
  }

  void ImgData::SurfTriangleSlot(bool val)
  {
    if(val and mesh->surfView() != Mesh::SURF_TRIANGLE) {
      emit changedInterface();
      mesh->setShowTriangle();
      emit viewerUpdate();
    }
  }

  void ImgData::SurfLabelSlot(bool val)
  {
    if(val and mesh->surfView() != Mesh::SURF_LABEL) {
      emit changedInterface();
      mesh->setShowLabel();
      emit viewerUpdate();
    }
  }

  void ImgData::VertexViewSlot(const QString &s)
  {
    if(s != mesh->surfVertexView()) {
      emit changedInterface();
      mesh->setShowVertex(s,false);
      emit viewerUpdate();
    }
  }

  void ImgData::TriangleViewSlot(const QString &s)
  {
    if(s != mesh->surfTriangleView()) {
      //emit changedInterface();
      mesh->setShowTriangle(s,false);
      updateTriColor();
      emit viewerUpdate();
    }
  }

  void ImgData::LabelViewSlot(const QString &s)
  {
    if(s != mesh->surfLabelView()) {
      bool needUpdateTris = false;
      emit changedInterface();
      // These require an update of the VBO
      if(s == "Wall Heat" or mesh->surfLabelView() == "Wall Heat"
         or s == "Cell Color" or mesh->surfLabelView() == "Cell Color")
        needUpdateTris = true;

      mesh->setShowLabel(s,false);
      // RSS Note this can not run in parallel with processes 
      // so this control needs to be disabled when processes are running
      if(needUpdateTris)
        updateTriColor();
      emit viewerUpdate();
    }
  } 

  void ImgData::VertexColorMapSlot()
  {
    if(mesh->surfVertexView() == "Signal")
      editSurfTransferFunction();
    else 
      return;

    emit viewerUpdate();
  } 

  void ImgData::TriangleColorMapSlot()
  {
    if(mesh->surfTriangleView() == "Triangle Value") {
      ColorMap triValueColors(&mesh->triangleValueColors(), &mesh->triangleValueBounds());
      triangleColorEditDlg = new ColorEditDlg(triValueColors, parent);
    } else if(mesh->surfTriangleView() == "Triangle Index") {
      ColorMap triIndexColors(&mesh->triangleIndexColors());
      triangleColorEditDlg = new ColorEditDlg(triIndexColors, parent);
    } else
      return;

    connect(triangleColorEditDlg, SIGNAL(update()), this, SLOT(UpdateColorMapSlot()));
    if(triangleColorEditDlg->exec() == QDialog::Accepted)
      emit viewerUpdate();
    delete triangleColorEditDlg;
  } 

  void ImgData::LabelColorMapSlot()
  {
    if(mesh->surfLabelView() == "Label Heat" or mesh->surfLabelView() == "Wall Heat")
      editHeatTransferFunction();
    else if(mesh->surfLabelView() == "Label")
      emit toggleEditLabels();
    else
      return;

    emit viewerUpdate();
  } 

  void ImgData::SurfParentSlot(bool val)
  {
    if(val != mesh->useParents()) {
      emit changedInterface();
      mesh->setUseParents(val);
      updateTriColor();
      emit viewerUpdate();
    }
  }
  
  void ImgData::SurfCullSlot(bool val)
  {
    if(val != mesh->culling()) {
      emit changedInterface();
      mesh->setCulling(val);
      emit viewerUpdate();
    }
  }
  
  // Mesh
  void ImgData::MeshShowSlot(bool val)
  {
    if(val != mesh->showMesh()) {
      emit changedInterface();
      mesh->setShowMesh(val);
      emit viewerUpdate();
    }
  }
   
  void ImgData::AxisShowSlot(bool val)
  {
    if(val != mesh->showAxis()) {
      emit changedInterface();
      mesh->setShowAxis(val);
      emit viewerUpdate();
    }
  } 

  void ImgData::ChangeMeshViewModeSlot(const QString &s)
  {
    if(s != mesh->meshView()) {
      mesh->setMeshView(s,false);
      emit viewerUpdate();
      emit changedInterface();
    }
  }

  void ImgData::ChangeAxisViewModeSlot(const QString &s)
  {
    if(s != mesh->axisView()) {
      mesh->setAxisView(s,false);
      emit viewerUpdate();
      emit changedInterface();
    }
  }
   
  void ImgData::MeshLinesSlot(bool val)
  {
    if(val != mesh->showMeshLines()) {
      emit changedInterface();
      mesh->setShowMeshLines(val);
      emit viewerUpdate();
    }
  }
  
  void ImgData::MeshPointsSlot(bool val)
  {
    if(val != mesh->showMeshPoints()) {
      emit changedInterface();
      mesh->setShowMeshPoints(val);
      emit viewerUpdate();
    }
  }
  
  void ImgData::CellMapSlot(bool val)
  {
    if(val != mesh->showMeshCellMap()) {
      emit changedInterface();
      mesh->setShowMeshCellMap(val);
      emit viewerUpdate();
    }
  }
  
  void ImgData::ShowTransSlot(bool val)
  {
    if(val != stack->showTrans()) {
      emit changedInterface();
      stack->setShowTrans(val);
      emit viewerUpdate();
    }
  }
  
  void ImgData::ShowBBoxSlot(bool val)
  {
    if(val != stack->showBBox()) {
      emit changedInterface();
      stack->setShowBBox(val);
      emit viewerUpdate();
    }
  }
  
  void ImgData::ShowScaleSlot(bool val)
  {
    if(val != stack->showScale()) {
      emit changedInterface();
      stack->setShowScale(val);
      updateStackSize();
      emit viewerUpdate();
    }
  }
  
  void ImgData::TieScalesSlot(bool val)
  {
    if(val != stack->tieScales()) {
      emit changedInterface();
      stack->setTieScales(val);
      emit viewerUpdate();
    }
  }

  void ImgData::ViewerUpdateSlot(void)
  {
    emit viewerUpdate();
  } 

  void ImgData::UpdateColorMapSlot(void)
  {
    updateTriColor();
    emit viewerUpdate();
  } 

  // Map back and forth between slider and scale
  int ImgData::toSliderScale(float s) 
  {
    return int((s - 1.0f) * 15000.0f);
  }
  
  float ImgData::fromSliderScale(int i) 
  {
    return 1.0f + float(i) / 15000.0f;
  }
  
  // Scale slots in x,y,z
  void ImgData::ScaleSlotX(int val)
  {
    float newScale = fromSliderScale(val);
    if(newScale != stack->scale().x()) {
      emit changedInterface();
      if(!stack->tieScales())
        stack->setScale(Point3f(newScale, stack->scale().y(), stack->scale().z()));
      else
        stack->setScale(Point3f(newScale, newScale, newScale));
      updateStackSize();
      emit updateSliderScale();
      emit viewerUpdate();
    }
  }

  void ImgData::ScaleSlotY(int val)
  {
    float newScale = fromSliderScale(val);
    if(newScale != stack->scale().y()) {
      emit changedInterface();
      if(!stack->tieScales())
        stack->setScale(Point3f(stack->scale().x(), newScale, stack->scale().z()));
      else
        stack->setScale(Point3f(newScale, newScale, newScale));
      updateStackSize();
      emit updateSliderScale();
      emit viewerUpdate();
    }
  }

  void ImgData::ScaleSlotZ(int val)
  {
    float newScale = fromSliderScale(val);
    if(newScale != stack->scale().z()) {
      emit changedInterface();
      if(!stack->tieScales())
        stack->setScale(Point3f(stack->scale().x(), stack->scale().y(), newScale));
      else
        stack->setScale(Point3f(newScale, newScale, newScale));
      updateStackSize();
      emit updateSliderScale();
      emit viewerUpdate();
    }
  }

  Point3f ImgData::imageGradientI(Point3i ipos)
  {
    uint ix = ipos.x();
    uint iy = ipos.y();
    uint iz = ipos.z();
    Point3f grad(0, 0, 0);
    const Point3f size(stack->size());
    double xScale = (size.x() / stack->size().x()) * 2.0;
    double yScale = (size.y() / stack->size().y()) * 2.0;
    double zScale = (size.z() / stack->size().z()) * 2.0;
    HVecUS& Data = currentData();
    grad[0] = (Data[offset(ix, iy, iz) + 1] - Data[offset(ix, iy, iz) - 1]) / xScale;
    grad[1] = (Data[offset(ix, iy, iz) + stack->size().x()] - Data[offset(ix, iy, iz) - stack->size().x()]) / yScale;
    grad[2] = (Data[offset(ix, iy, iz) + stack->size().x() * stack->size().y()]
               - Data[offset(ix, iy, iz) - stack->size().x() * stack->size().y()]) / zScale;
    return grad;
  }
  
  Point3f ImgData::imageGradientW(Point3d worldpos)
  {
    Point3i ipos(worldToImagei(Point3f(worldpos)));
    return imageGradientI(ipos);
  }
  
  uint ImgData::imageLevel(Point3d worldpos)
  {
    HVecUS& Data = currentData();
    Point3i ipos = worldToImagei(Point3f(worldpos));
    return Data[offset(ipos)];
  }
  
  Point2i ImgData::imageMinMax()
  {
    HVecUS& Data = currentData();
    uint min = 65535;
    uint max = 0;
    for(size_t i = 0; i < stack->storeSize(); ++i) {
      uint val = Data[i];
      if(val > max)
        max = val;
      if(val < min)
        min = val;
    }
    return Point2i(min, max);
  }
}
