//
// 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 "MorphoViewer.hpp"

#include "cuda/CudaExport.hpp"
#include "CutSurf.hpp"
#include "Dir.hpp"
#include "ImageData.hpp"
#include "Information.hpp"
#include "MorphoGraphX.hpp"
#include "Process.hpp"

#include <iostream>
#include <QDomElement>
#include <QFileInfo>
#include <QGLFormat>
#include <QGLFramebufferObject>
#include <QGLFramebufferObjectFormat>
#include <QImage>
#include <QShowEvent>
#include <QMessageBox>
#include <QFileDialog>
#include <string>

#ifdef EPSILON
#  undef EPSILON
#endif
#define EPSILON 1e-3

using namespace mgx;
using namespace qglviewer;

MGXKeyFrameInterpolator::MGXKeyFrameInterpolator(Frame* fr) : KeyFrameInterpolator(fr) {}

void MGXKeyFrameInterpolator::interpolateAtTime(float time) 
{
  KeyFrameInterpolator::interpolateAtTime(time);
}

MGXCameraFrame::MGXCameraFrame() : qglviewer::ManipulatedCameraFrame(), _zoom(0) {}

void MGXCameraFrame::wheelEvent(QWheelEvent* const event, Camera* const camera)
{
  const float wheelSensitivityCoef = 8E-4f;
  switch(action_) {
  case QGLViewer::ZOOM:
    _zoom -= wheelSensitivity() * event->delta() * wheelSensitivityCoef;
    if(DEBUG)
      Information::out << "  zoom = " << _zoom << endl;
    emit manipulated();
    break;
  default:
    wheelEvent(event, camera);
    return;
  }

  if(previousConstraint_)
    setConstraint(previousConstraint_);

  const int finalDrawAfterWheelEventDelay = 400;

  QTimer::singleShot(finalDrawAfterWheelEventDelay, this, SLOT(flyUpdate()));

  // This could also be done *before* manipulated is emitted, so that isManipulated() returns false.
  // But then fastDraw would not be used with wheel.
  // Detecting the last wheel event and forcing a final draw() is done using the timer_.
  action_ = QGLViewer::NO_MOUSE_ACTION;
}

MGXCamera::MGXCamera() : qglviewer::Camera()
{
  setType(qglviewer::Camera::ORTHOGRAPHIC);
  _frame = new MGXCameraFrame();
  setFrame(_frame);
}

void MGXCamera::getOrthoWidthHeight(GLdouble& halfWidth, GLdouble& halfHeight) const
{
  qglviewer::Camera::getOrthoWidthHeight(halfWidth, halfHeight);
  float z = 0.5f * exp(_frame->zoom());
  halfWidth *= z;
  halfHeight *= z;
}

void MGXCamera::resetZoom() { _frame->setZoom(0.0f); }

void MGXCamera::fitSphere(const Vec& center, float radius)
{
  _frame->setZoom(0.0f);
  qglviewer::Camera::fitSphere(center, radius);
}

void MGXCamera::addKeyFrameToPath(int i) 
{
  Camera::addKeyFrameToPath(i);
}

void MGXCamera::playPath(int i) 
{
  Camera::playPath(i);
}

void MGXCamera::deletePath(int i) 
{
  Camera::deletePath(i);
}

void MGXCamera::resetPath(int i) 
{
  Camera::resetPath(i);
}

void MGXCamera::drawAllPaths() 
{
  Camera::drawAllPaths();
}

MorphoViewer::MorphoViewer(QWidget* parent)
  : QGLViewer(parent), quitting(false), DrawCellMap(false), DrawClipBox(false), 
    voxelEditCursor(false), pixelRadius(1), shiftPressed(false), altPressed(false), 
    controlPressed(false), leftButton(false), rightButton(false), guiAction(MESH_PICK_LABEL), 
    guiActionOn(false), meshPickFill(true), stackPickFill(true), stack1(NULL), stack2(NULL), cutSurf(NULL), 
    FlySpeed(0.005f), selectedLabel(0), sampling(2), fast_draw(false), MaxNbPeels(20), 
    GlobalBrightness(0.0), GlobalContrast(1.0), UnsharpStrength(1.0), Spinning(10000), 
    _camera(new MGXCamera()), SceneRadius(1e7), initialized(false), lastInitCamera(0),
    _processRunning(false)
{
  initObject(parent);
  setCamera(_camera);
}

void MorphoViewer::initObject(QWidget*)
{
  show_slice = -1;

  setAnimationPeriod(0);

  // Remap movement to manipulated frame
  setMouseBinding(Qt::SHIFT + Qt::LeftButton, CAMERA, TRANSLATE);
  setMouseBinding(Qt::CTRL + Qt::SHIFT + Qt::LeftButton, FRAME, TRANSLATE);
  // Prevent manipulated frame zooming
  setWheelBinding(Qt::NoModifier, CAMERA, ZOOM);
  // Set that to move along the z axis
  setMouseBinding(Qt::CTRL + Qt::MidButton, FRAME, MOVE_FORWARD);

  // Include shift in selection
  setShortcut(SAVE_SCREENSHOT, Qt::SHIFT + Qt::Key_Print);
  setShortcut(CAMERA_MODE, 0);
  setShortcut(EXIT_VIEWER, 0);
  setShortcut(INCREASE_FLYSPEED, 0);
  setShortcut(DECREASE_FLYSPEED, 0);
  setContextMenuPolicy(Qt::ActionsContextMenu);
  QGLFormat f = format();

  // Enable capture of mouse move events
  setMouseTracking(true);

  // Hotkey camera locations

  // Turn off buffer swap, we'll do it ourselves
  setAutoBufferSwap(false);
}

MorphoViewer::~MorphoViewer()
{
  saveStateToFile();

  // Clean FBO

  glDeleteTextures(NB_FRAMES, colorTexId);
  glDeleteTextures(NB_FRAMES, depthTexId);
  for(int i = 0; i < NB_FRAMES; ++i) {
    colorTexId[i] = 0;
    depthTexId[i] = 0;
  }
  glfuncs->glDeleteRenderbuffers(1, &depthBuffer);
  depthBuffer = 0;

  if(fboId != 0)
    glfuncs->glDeleteFramebuffers(1, &fboId);
  fboId = 0;
  if(fboCopyId != 0)
    glfuncs->glDeleteFramebuffers(1, &fboCopyId);
  fboCopyId = 0;
}

void MorphoViewer::init()
{
  // Check multi-sampling is de-activate, but in practice it seems useless ..
  if(format().sampleBuffers()) {
    QMessageBox::critical(this, 
      "There is a Problem with the OpenGL context", 
      "Multi-sampling is activated in the OpenGL context.\n"
      "Make sure your drivers are not setup to override application preferences.\n"
      "Note that as long as this persists, selection will not work properly in MorphoGraphX.");
  }

  QGLViewer::init();

  glfuncs = this->context()->contextHandle()->versionFunctions<MGXOpenGLFunctions>();
  if(!glfuncs ||
     !(glfuncs->initializeOpenGLFunctions())) {
    Information::out << "Error initializing OpenGL functions, MorphoGraphX requires OpenGL "
                     << MGX_REQUIRED_OPENGL_VERSION << endl;
    QCoreApplication::exit(126);
    return;
  }
  Information::out << "Rendering OpenGL " << (char*)(glGetString(GL_VERSION))
                   << " on " << (char*)(glGetString(GL_RENDERER)) << endl;

  GLint maxTexSize;
  glfuncs->glGetIntegerv(GL_MAX_TEXTURE_SIZE, &maxTexSize);
  Information::out << "OpenGL maximum texture size: " << maxTexSize << endl;

  initCamera();

  raycastingShader1.setVerbosity(0);
  raycastingShader1.init();

  raycastingShader1.addVertexShader(":/shaders/RayCasting.vert");
  raycastingShader1.addFragmentShader(":/shaders/Utils.frag");
  raycastingShader1.addFragmentShader(":/shaders/3DColoring.frag");
  raycastingShader1.addFragmentShaderCode("");   // To be replaced by the needed function
  raycastingShader1.addFragmentShader(":/shaders/RayCasting.frag");

  raycastingShader2.setVerbosity(0);
  raycastingShader2.init();

  raycastingShader2.addVertexShader(":/shaders/RayCasting.vert");
  raycastingShader2.addFragmentShader(":/shaders/Utils.frag");
  raycastingShader2.addFragmentShader(":/shaders/3DColoring.frag");
  raycastingShader2.addFragmentShaderCode("");   // To be replaced by the needed function
  raycastingShader2.addFragmentShader(":/shaders/RayCasting.frag");

  finalCombineShader.setVerbosity(0);
  finalCombineShader.init();

  finalCombineShader.addVertexShader(":/shaders/Combine.vert");
  finalCombineShader.addFragmentShader(":/shaders/Utils.frag");
  finalCombineShader.addFragmentShader(":/shaders/FinalCombine.frag");

  combineShader.setVerbosity(0);
  combineShader.init();

  combineShader.addVertexShader(":/shaders/Combine.vert");
  combineShader.addFragmentShader(":/shaders/Utils.frag");
  combineShader.addFragmentShader(":/shaders/Combine.frag");

  occlusionShader.setVerbosity(0);
  occlusionShader.init();

  occlusionShader.addVertexShader(":/shaders/PostProcess.vert");
  occlusionShader.addFragmentShader(":/shaders/Occlusion.frag");

  postProcessShader.setVerbosity(0);
  postProcessShader.init();

  postProcessShader.addVertexShader(":/shaders/PostProcess.vert");
  postProcessShader.addFragmentShader(":/shaders/PostProcess.frag");

  renderDepthShader.setVerbosity(0);
  renderDepthShader.init();

  renderDepthShader.addVertexShader(":/shaders/RenderDepth.vert");
  renderDepthShader.addFragmentShader(":/shaders/RenderDepth.frag");

  textureSurfShader.setVerbosity(0);
  textureSurfShader.init();

  textureSurfShader.addVertexShader(":/shaders/Light.vert");
  textureSurfShader.addVertexShader(":/shaders/TextureSurf.vert");
  textureSurfShader.addVertexShader(":/shaders/PeelingSurf.vert");
  textureSurfShader.addFragmentShader(":/shaders/Light.frag");
  textureSurfShader.addFragmentShader(":/shaders/TextureSurf.frag");
  textureSurfShader.addFragmentShader(":/shaders/PeelingSurf.frag");

  volumeSurfShader1.setVerbosity(0);
  volumeSurfShader1.init();

  volumeSurfShader1.addVertexShader(":/shaders/Light.vert");
  volumeSurfShader1.addVertexShader(":/shaders/VolumeSurf.vert");
  volumeSurfShader1.addVertexShader(":/shaders/PeelingSurf.vert");
  volumeSurfShader1.addFragmentShader(":/shaders/Light.frag");
  volumeSurfShader1.addFragmentShader(":/shaders/Utils.frag");
  volumeSurfShader1.addFragmentShader(":/shaders/3DColoring.frag");
  volumeSurfShader1.addFragmentShaderCode("");
  volumeSurfShader1.addFragmentShader(":/shaders/VolumeSurf.frag");
  volumeSurfShader1.addFragmentShader(":/shaders/PeelingSurf.frag");

  volumeSurfShader2.setVerbosity(0);
  volumeSurfShader2.init();

  volumeSurfShader2.addVertexShader(":/shaders/Light.vert");
  volumeSurfShader2.addVertexShader(":/shaders/VolumeSurf.vert");
  volumeSurfShader2.addVertexShader(":/shaders/PeelingSurf.vert");
  volumeSurfShader2.addFragmentShader(":/shaders/Light.frag");
  volumeSurfShader2.addFragmentShader(":/shaders/Utils.frag");
  volumeSurfShader2.addFragmentShader(":/shaders/3DColoring.frag");
  volumeSurfShader2.addFragmentShaderCode("");
  volumeSurfShader2.addFragmentShader(":/shaders/VolumeSurf.frag");
  volumeSurfShader2.addFragmentShader(":/shaders/PeelingSurf.frag");

  indexSurfShader.setVerbosity(0);
  indexSurfShader.init();

  indexSurfShader.addVertexShader(":/shaders/Light.vert");
  indexSurfShader.addVertexShader(":/shaders/IndexSurf.vert");
  indexSurfShader.addVertexShader(":/shaders/PeelingSurf.vert");
  indexSurfShader.addFragmentShader(":/shaders/Light.frag");
  indexSurfShader.addFragmentShader(":/shaders/IndexSurf.frag");
  indexSurfShader.addFragmentShader(":/shaders/PeelingSurf.frag");

  colorSurfShader.setVerbosity(0);
  colorSurfShader.init();

  colorSurfShader.addVertexShader(":/shaders/Light.vert");
  colorSurfShader.addVertexShader(":/shaders/ColorSurf.vert");
  colorSurfShader.addVertexShader(":/shaders/PeelingSurf.vert");
  colorSurfShader.addFragmentShader(":/shaders/Light.frag");
  colorSurfShader.addFragmentShader(":/shaders/ColorSurf.frag");
  colorSurfShader.addFragmentShader(":/shaders/PeelingSurf.frag");

  REPORT_GL_ERROR("OpenGL error initializing shaders");

  // camera()->setType(qglviewer::Camera::ORTHOGRAPHIC);
  camera()->setSceneRadius(SceneRadius);
  camera()->setSceneCenter(Vec(0, 0, 0));
  camera()->centerScene();
  camera()->showEntireScene();

  // Color3f clearColor = Colors::getColor(Colors::BackgroundColor);
  // glClearColor(clearColor.r(), clearColor.g(), clearColor.b(), 1.0);
  // glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT | GL_ACCUM_BUFFER_BIT);
  glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, GL_TRUE);
  glShadeModel(GL_SMOOTH);
  // glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
  // glTexEnvi(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_BLEND);
  glEnable(GL_DEPTH_CLAMP_NV);   // Allow to zoom more, but cause reverse later. DO NOT REMOVE

  stack1->initTex();
  stack2->initTex();
  ImgData::scaleBar.init(this);

  // Init FBO
  glfuncs->glGenFramebuffers(1, &fboId);
  glfuncs->glGenFramebuffers(1, &fboCopyId);

  // Create textures for depth and color buffer
  depthBuffer = 0;

  for(int i = 0; i < NB_FRAMES; ++i) {
    colorTexId[i] = 0;
    depthTexId[i] = 0;
  }

  prevWidth = prevHeight = 0;
  // updateFBOTex(width(), height());

  // To start with, we always render on the screen
  baseFboId = 0;
  //  Information::out << "OpenGL initialized" << endl;
}

void MorphoViewer::postDraw()
{
  if(quitting)
    return;
  glEnable(GL_DEPTH_TEST);
  QGLViewer::postDraw();
  swapBuffers();
}

void MorphoViewer::initFromDOMElement(const QDomElement& element)
{
  // Restore standard state
  QGLViewer::initFromDOMElement(element);

  QGLFormat def = QGLFormat::defaultFormat();

  QDomElement child = element.firstChild().toElement();
  while(!child.isNull()) {
    if(child.tagName() == "GLFormat") {
      if(child.hasAttribute("alpha"))
        def.setAlpha(child.attribute("alpha").toLower() == "yes");
      if(child.hasAttribute("depth"))
        def.setDepth(child.attribute("depth").toLower() == "yes");
      if(child.hasAttribute("rgba"))
        def.setRgba(child.attribute("rgba").toLower() == "yes");
      if(child.hasAttribute("stereo"))
        def.setStereo(child.attribute("stereo").toLower() == "yes");
      if(child.hasAttribute("stencil"))
        def.setStencil(child.attribute("stencil").toLower() == "yes");
      if(child.hasAttribute("doubleBuffer"))
        def.setDoubleBuffer(child.attribute("doubleBuffer").toLower() == "yes");
      if(child.hasAttribute("sampleBuffers"))
        def.setSampleBuffers(child.attribute("sampleBuffers").toLower() == "yes");
      if(child.hasAttribute("directRendering"))
        def.setDirectRendering(child.attribute("directRendering").toLower() == "yes");
      if(child.hasAttribute("hasOverlay"))
        def.setOverlay(child.attribute("hasOverlay").toLower() == "yes");
      if(def.accum())
        if(child.hasAttribute("accumBufferSize"))
          def.setAccumBufferSize(child.attribute("accumBufferSize").toInt());
      if(def.alpha())
        if(child.hasAttribute("alphaBufferSize"))
          def.setAlphaBufferSize(child.attribute("alphaBufferSize").toInt());
      if(def.depth())
        if(child.hasAttribute("depthBufferSize"))
          def.setDepthBufferSize(child.attribute("depthBufferSize").toInt());
      if(def.sampleBuffers())
        if(child.hasAttribute("samples"))
          def.setSamples(child.attribute("samples").toInt());
      if(def.stencil())
        if(child.hasAttribute("stencilBufferSize"))
          def.setStencilBufferSize(child.attribute("stencilBufferSize").toInt());
      if(def.rgba()) {
        if(child.hasAttribute("redBufferSize"))
          def.setRedBufferSize(child.attribute("redBufferSize").toInt());
        if(child.hasAttribute("greenBufferSize"))
          def.setGreenBufferSize(child.attribute("greenBufferSize").toInt());
        if(child.hasAttribute("blueBufferSize"))
          def.setBlueBufferSize(child.attribute("blueBufferSize").toInt());
      }
    }
    child = child.nextSibling().toElement();
  }
  QGLFormat::setDefaultFormat(def);
}

QDomElement MorphoViewer::domElement(const QString& name, QDomDocument& document) const
{
  // Creates a custom node for a light
  QDomElement de = document.createElement("GLFormat");
  QGLFormat def = QGLFormat::defaultFormat();
  de.setAttribute("accum", (def.accum() ? "yes" : "no"));
  if(def.accum() && def.accumBufferSize() != -1)
    de.setAttribute("accumBufferSize", std::max(0, def.accumBufferSize()));
  de.setAttribute("alpha", (def.alpha() ? "yes" : "no"));
  if(def.alpha() && def.alphaBufferSize() != -1)
    de.setAttribute("alphaBufferSize", std::max(0, def.alphaBufferSize()));
  de.setAttribute("depth", (def.depth() ? "yes" : "no"));
  if(def.depth() && def.depthBufferSize() != -1)
    de.setAttribute("depthBufferSize", std::max(0, def.depthBufferSize()));
  de.setAttribute("doubleBuffer", (def.doubleBuffer() ? "yes" : "no"));
  de.setAttribute("directRendering", (def.directRendering() ? "yes" : "no"));
  de.setAttribute("hasOverlay", (def.hasOverlay() ? "yes" : "no"));
  de.setAttribute("rgba", (def.rgba() ? "yes" : "no"));
  if(def.rgba()) {
    if(def.redBufferSize() != -1)
      de.setAttribute("redBufferSize", std::max(0, def.redBufferSize()));
    if(def.greenBufferSize() != -1)
      de.setAttribute("greenBufferSize", std::max(0, def.greenBufferSize()));
    if(def.blueBufferSize() != -1)
      de.setAttribute("blueBufferSize", std::max(0, def.blueBufferSize()));
  }
  de.setAttribute("sampleBuffers", (def.sampleBuffers() ? "yes" : "no"));
  if(def.sampleBuffers() && def.samples() != -1)
    de.setAttribute("samples", std::max(0, def.samples()));
  de.setAttribute("stereo", (def.stereo() ? "yes" : "no"));
  de.setAttribute("stencil", (def.stencil() ? "yes" : "no"));
  if(def.stencil() && def.stencilBufferSize() != -1)
    de.setAttribute("stencilBufferSize", std::max(0, def.stencilBufferSize()));

  // Get default state domElement and append custom node
  QDomElement res = QGLViewer::domElement(name, document);
  res.appendChild(de);
  return res;
}

void MorphoViewer::initFormat()
{
  QGLFormat def = QGLFormat::defaultFormat();
  QFile file(".qglviewer.xml");
  if(file.open(QIODevice::ReadOnly)) {
    QDomDocument doc("QGLViewer");
    doc.setContent(&file);
    file.close();
    QDomElement root = doc.firstChildElement("QGLViewer");
    QDomElement child = root.firstChildElement("GLFormat");
    if(!child.isNull()) {
      if(child.hasAttribute("alpha"))
        def.setAlpha(child.attribute("alpha").toLower() == "yes");
      if(child.hasAttribute("depth"))
        def.setDepth(child.attribute("depth").toLower() == "yes");
      if(child.hasAttribute("rgba"))
        def.setRgba(child.attribute("rgba").toLower() == "yes");
      if(child.hasAttribute("stereo"))
        def.setStereo(child.attribute("stereo").toLower() == "yes");
      if(child.hasAttribute("stencil"))
        def.setStencil(child.attribute("stencil").toLower() == "yes");
      if(child.hasAttribute("doubleBuffer"))
        def.setDoubleBuffer(child.attribute("doubleBuffer").toLower() == "yes");
      if(child.hasAttribute("sampleBuffers"))
        def.setSampleBuffers(child.attribute("sampleBuffers").toLower() == "yes");
      if(child.hasAttribute("directRendering"))
        def.setDirectRendering(child.attribute("directRendering").toLower() == "yes");
      if(child.hasAttribute("hasOverlay"))
        def.setOverlay(child.attribute("hasOverlay").toLower() == "yes");
      if(def.accum())
        if(child.hasAttribute("accumBufferSize"))
          def.setAccumBufferSize(child.attribute("accumBufferSize").toInt());
      if(def.alpha())
        if(child.hasAttribute("alphaBufferSize"))
          def.setAlphaBufferSize(child.attribute("alphaBufferSize").toInt());
      if(def.depth())
        if(child.hasAttribute("depthBufferSize"))
          def.setDepthBufferSize(child.attribute("depthBufferSize").toInt());
      if(def.sampleBuffers())
        if(child.hasAttribute("samples"))
          def.setSamples(child.attribute("samples").toInt());
      if(def.stencil())
        if(child.hasAttribute("stencilBufferSize"))
          def.setStencilBufferSize(child.attribute("stencilBufferSize").toInt());
      if(def.rgba()) {
        if(child.hasAttribute("redBufferSize"))
          def.setRedBufferSize(child.attribute("redBufferSize").toInt());
        if(child.hasAttribute("greenBufferSize"))
          def.setGreenBufferSize(child.attribute("greenBufferSize").toInt());
        if(child.hasAttribute("blueBufferSize"))
          def.setBlueBufferSize(child.attribute("blueBufferSize").toInt());
      }
    }
  }
  QGLFormat::setDefaultFormat(def);
}

void MorphoViewer::showEvent(QShowEvent* event)
{
  return;
  QGLViewer::showEvent(event);
  if(!initialized) {
    restoreStateFromFile();
    lastInitCamera = camera();
    initialized = true;
  }
}

void MorphoViewer::initCamera()
{
  if(camera() != lastInitCamera) {
    restoreStateFromFile();
    lastInitCamera = camera();
  }
}

void MorphoViewer::preDraw()
{
  if(quitting)
    return;
  QGLViewer::preDraw();
  camera()->setSceneRadius(SceneRadius);

  glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST);

  if(!finalCombineShader.initialized())
    if(!finalCombineShader.setupShaders()) {
      Information::err << "Error compiling combine shaders" << endl;
    }
  if(!combineShader.initialized())
    if(!combineShader.setupShaders()) {
      Information::err << "Error compiling combine shaders" << endl;
    }
  if(!occlusionShader.initialized())
    if(!occlusionShader.setupShaders()) {
      Information::err << "Error compiling combine shaders" << endl;
    }
  if(!postProcessShader.initialized())
    if(!postProcessShader.setupShaders()) {
      Information::err << "Error compiling combine shaders" << endl;
    }
  if(!renderDepthShader.initialized())
    if(!renderDepthShader.setupShaders()) {
      Information::err << "Error compiling render_depth shaders" << endl;
    }
  stack1->setupVolumeShader(raycastingShader1, 2);
  if(!raycastingShader1.initialized())
    if(!raycastingShader1.setupShaders()) {
      Information::err << "Error compiling colormap shaders" << endl;
    }
  stack2->setupVolumeShader(raycastingShader2, 2);
  if(!raycastingShader2.initialized())
    if(!raycastingShader2.setupShaders()) {
      Information::err << "Error compiling colormap shaders" << endl;
    }
  if(!textureSurfShader.initialized())
    if(!textureSurfShader.setupShaders()) {
      Information::err << "Error compiling texture_surf shaders" << endl;
    }
  stack1->setupVolumeShader(volumeSurfShader1, 3);
  // if(DEBUG)
  // Information::out << "code for volumeSurfShader1 color:\n" <<
  // volumeSurfShader1.getFragmentShader(2).toStdString() << endl;
  if(!volumeSurfShader1.initialized())
    if(!volumeSurfShader1.setupShaders()) {
      Information::err << "Error compiling volume_surf shaders" << endl;
    }
  stack2->setupVolumeShader(volumeSurfShader2, 3);
  // if(DEBUG)
  // Information::out << "code for volumeSurfShader2 color:\n" <<
  // volumeSurfShader2.getFragmentShader(2).toStdString() << endl;
  if(!volumeSurfShader2.initialized())
    if(!volumeSurfShader2.setupShaders()) {
      Information::err << "Error compiling volume_surf shaders" << endl;
    }
  if(!indexSurfShader.initialized())
    if(!indexSurfShader.setupShaders()) {
      Information::err << "Error compiling index_surf shaders" << endl;
    }
  if(!colorSurfShader.initialized())
    if(!colorSurfShader.setupShaders()) {
      Information::err << "Error compiling color_surf shaders" << endl;
    }
}

static bool checkFBO(const char* file, size_t ln)
{
  GLenum status = glfuncs->glCheckFramebufferStatus(GL_FRAMEBUFFER_EXT);

  if(status == 0) {
    Information::err << "Error line " << ln << " creating framebuffer" << endl;
    return false;
  } else {
    if(status != GL_FRAMEBUFFER_COMPLETE_EXT) {
      Information::err << "Error file " << file << " on line " << ln << ": ";
      switch(status) {
      case GL_FRAMEBUFFER_INCOMPLETE_ATTACHMENT_EXT:
        Information::err << "Framebuffer is incomplete: incomplete attachment" << endl;
        break;
      case GL_FRAMEBUFFER_INCOMPLETE_MISSING_ATTACHMENT_EXT:
        Information::err << "Framebuffer is incomplete: missing attachment" << endl;
        break;
      case GL_FRAMEBUFFER_INCOMPLETE_DRAW_BUFFER_EXT:
        Information::err << "Framebuffer is incomplete: incomplete draw buffer" << endl;
        break;
      case GL_FRAMEBUFFER_INCOMPLETE_READ_BUFFER_EXT:
        Information::err << "Framebuffer is incomplete: incomplete read buffer" << endl;
        break;
      case GL_FRAMEBUFFER_INCOMPLETE_FORMATS_EXT:
        Information::err << "Framebuffer is incomplete: incomplete formats" << endl;
        break;
      case GL_FRAMEBUFFER_INCOMPLETE_DIMENSIONS_EXT:
        Information::err << "Framebuffer is incomplete: incomplete dimensions" << endl;
        break;
      case GL_FRAMEBUFFER_UNSUPPORTED_EXT:
        Information::err << "Framebuffer unsupported" << endl;
        break;
      default:
        Information::err << "Framebuffer is incomplete: status = " << status << endl;
      }
      return false;
    }
  }
  return true;
}

void MorphoViewer::ScreenSamplingSlot(int val) {
  sampling = val / 10.0f + 1.0f;
}

void MorphoViewer::updateFBOTex(int width, int height, bool fast_draw)
{
  if(DEBUG)
    Information::out << "updateFBOTex(" << width << ", " << height << ", " << fast_draw << ")" << endl;
  if(fast_draw) {
    texWidth = int(ceil(width / sampling));
    texHeight = int(ceil(height / sampling));
  } else {
    texWidth = width;
    texHeight = height;
  }
  if(texWidth == prevWidth and texHeight == prevHeight)
    return;
  // Information::err << "width = " << texWidth << " - height = " << texHeight << endl;

  prevWidth = texWidth;
  prevHeight = texHeight;

  if(colorTexId[0] == 0)
    glGenTextures(NB_FRAMES, colorTexId);
  if(depthTexId[0] == 0)
    glGenTextures(NB_FRAMES, depthTexId);
  if(depthBuffer == 0)
    glfuncs->glGenRenderbuffers(1, &depthBuffer);

  // Create textures
  for(int i = 0; i < NB_FRAMES; ++i) {
    glBindTexture(GL_TEXTURE_2D, colorTexId[i]);

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

    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, texWidth, texHeight, 0, GL_RGBA, GL_UNSIGNED_SHORT, NULL);

    glBindTexture(GL_TEXTURE_2D, depthTexId[i]);

    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
    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_COMPARE_MODE, GL_COMPARE_R_TO_TEXTURE);
    glTexParameteri(GL_TEXTURE_2D, GL_DEPTH_TEXTURE_MODE, GL_ALPHA);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_COMPARE_FUNC, GL_LEQUAL);

    glTexImage2D(GL_TEXTURE_2D, 0, GL_DEPTH_COMPONENT32F, texWidth, texHeight, 0, GL_DEPTH_COMPONENT, GL_FLOAT,
                 NULL);
  }

  glBindTexture(GL_TEXTURE_2D, 0);

  // Bind textures to framebuffer

  glfuncs->glBindFramebuffer(GL_FRAMEBUFFER_EXT, fboId);

  glfuncs->glBindRenderbuffer(GL_RENDERBUFFER_EXT, depthBuffer);
  glfuncs->glRenderbufferStorage(GL_RENDERBUFFER_EXT, GL_DEPTH_COMPONENT32F, texWidth, texHeight);
  glfuncs->glBindRenderbuffer(GL_RENDERBUFFER_EXT, 0);

  glfuncs->glFramebufferTexture2D(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, GL_TEXTURE_2D, colorTexId[0], 0);
  glfuncs->glFramebufferTexture2D(GL_FRAMEBUFFER_EXT, GL_DEPTH_ATTACHMENT_EXT, GL_TEXTURE_2D, depthTexId[0], 0);
  REPORT_GL_ERROR("CREATE_FBO");

  checkFBO(__FILE__, __LINE__);
}

void MorphoViewer::resizeGL(int width, int height)
{
  drawWidth = width;
  drawHeight = height;
  QGLViewer::resizeGL(width, height);
}

void MorphoViewer::setupCopyFB(GLuint depth, GLint color)
{
  glfuncs->glBindFramebuffer(GL_FRAMEBUFFER_EXT, fboCopyId);
  if(color != 0)
    glfuncs->glFramebufferTexture2D(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, GL_TEXTURE_2D, color, 0);
  if(depth != 0)
    glfuncs->glFramebufferTexture2D(GL_FRAMEBUFFER_EXT, GL_DEPTH_ATTACHMENT_EXT, GL_TEXTURE_2D, depth, 0);
}

void MorphoViewer::setupFramebuffer(GLuint depth, GLuint color, GLbitfield clear)
{
  GLdouble proj[16], mv[16];
  glGetDoublev(GL_MODELVIEW_MATRIX, mv);
  glGetDoublev(GL_PROJECTION_MATRIX, proj);

  /*
   * if(DEBUG)
   *{
   *  Information::out << "GL_PROJECTION_MATRIX = " << endl;
   *  for(int i = 0 ; i < 4 ; ++i)
   *  {
   *    Information::out << "{ ";
   *    for(int j = 0 ; j < 4 ; ++j)
   *    {
   *      Information::out << proj[i+4*j] << " ";
   *    }
   *    Information::out << "}\n";
   *  }
   *  Information::out << endl;
   *}
   */

  glfuncs->glBindFramebuffer(GL_FRAMEBUFFER_EXT, fboId);
  glfuncs->glFramebufferTexture2D(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, GL_TEXTURE_2D, color, 0);
  if(depth != 0) {
    glfuncs->glFramebufferTexture2D(GL_FRAMEBUFFER_EXT, GL_DEPTH_ATTACHMENT_EXT, GL_TEXTURE_2D, depth, 0);
  } else {
    glfuncs->glFramebufferRenderbuffer(GL_FRAMEBUFFER_EXT, GL_DEPTH_ATTACHMENT_EXT, GL_RENDERBUFFER_EXT, depthBuffer);
  }

  // checkFBO(__LINE__);

  glPushAttrib(GL_VIEWPORT_BIT);
  glViewport(0, 0, texWidth, texHeight);

  glClear(clear);

  glMatrixMode(GL_PROJECTION);
  glPushMatrix();
  glLoadIdentity();
  glMultMatrixd(proj);

  glMatrixMode(GL_MODELVIEW);
  glPushMatrix();
  glLoadIdentity();
  glMultMatrixd(mv);
}

void MorphoViewer::resetupFramebuffer(GLuint depth, GLuint color, GLbitfield clear)
{
  glfuncs->glFramebufferTexture2D(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, GL_TEXTURE_2D, color, 0);
  if(depth > 0) {
    glBindTexture(GL_TEXTURE_2D, depth);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_COMPARE_MODE, GL_COMPARE_R_TO_TEXTURE);
    glBindTexture(GL_TEXTURE_2D, 0);
    glfuncs->glFramebufferTexture2D(GL_FRAMEBUFFER_EXT, GL_DEPTH_ATTACHMENT_EXT, GL_TEXTURE_2D, depth, 0);
  } else {
    glfuncs->glFramebufferRenderbuffer(GL_FRAMEBUFFER_EXT, GL_DEPTH_ATTACHMENT_EXT, GL_RENDERBUFFER_EXT, depthBuffer);
  }

  glClear(clear);
}

void MorphoViewer::resetFramebuffer()
{
  glMatrixMode(GL_PROJECTION);
  glPopMatrix();
  glMatrixMode(GL_MODELVIEW);
  glPopMatrix();

  glPopAttrib();

  glfuncs->glBindFramebuffer(GL_FRAMEBUFFER_EXT, baseFboId);
}

void MorphoViewer::startScreenCoordinatesSystem(bool upward) const
{
  glMatrixMode(GL_PROJECTION);
  glPushMatrix();
  glLoadIdentity();
  if(upward)
    glOrtho(0, current_device->width(), 0, current_device->height(), 0.0, -1.0);
  else
    glOrtho(0, current_device->width(), current_device->height(), 0, 0.0, -1.0);

  glMatrixMode(GL_MODELVIEW);
  glPushMatrix();
  glLoadIdentity();
}

void MorphoViewer::drawColorTexture(int i, bool draw_depth)
{
  // Check the texture ...
  glPolygonMode(GL_FRONT, GL_FILL);
  glEnable(GL_TEXTURE_2D);
  glMatrixMode(GL_TEXTURE);
  glPushMatrix();
  glLoadIdentity();
  glMatrixMode(GL_MODELVIEW);

  if(i >= 0) {
    if(draw_depth)
      glBindTexture(GL_TEXTURE_2D, depthTexId[i]);
    else
      glBindTexture(GL_TEXTURE_2D, colorTexId[i]);
  }

  // For now, render the depth buffer instead of the color one
  startScreenCoordinatesSystem(true);

  glBegin(GL_QUADS);
  glColor4f(1, 1, 1, 1);
  glTexCoord2d(0, 0);
  glVertex2d(0, 0);
  glTexCoord2d(1, 0);
  glVertex2d(drawWidth, 0);
  glTexCoord2d(1, 1);
  glVertex2d(drawWidth, drawHeight);
  glTexCoord2d(0, 1);
  glVertex2d(0, drawHeight);
  glEnd();

  stopScreenCoordinatesSystem();

  glMatrixMode(GL_TEXTURE);
  glPopMatrix();
  glMatrixMode(GL_MODELVIEW);

  if(i >= 0)
    glBindTexture(GL_TEXTURE_2D, 0);
  glDisable(GL_TEXTURE_2D);
}

void MorphoViewer::alternatePeels(int& curPeelId, int& prevPeelId, int fullImgId)
{
  if(curPeelId == -1) {
    curPeelId = 3;
  } else if(curPeelId == 3) {
    prevPeelId = 3;
    curPeelId = 4;
  } else {
    prevPeelId = 4;
    curPeelId = 3;
  }

  if(fullImgId != -1) {
    Shader::activeTexture(Shader::AT_FRONT_TEX);
    glBindTexture(GL_TEXTURE_2D, depthTexId[fullImgId]);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_COMPARE_MODE, GL_NONE);
    Shader::activeTexture(Shader::AT_NONE);
  }

  resetupFramebuffer(depthTexId[curPeelId], colorTexId[curPeelId]);
}

void MorphoViewer::combinePeels(int& fullImgId, int curPeelId, int volume1, int volume2)
{
  int newFinalId = (fullImgId == FI_FULL_IMG1) ? FI_FULL_IMG2 : FI_FULL_IMG1;
  resetupFramebuffer(depthTexId[newFinalId], colorTexId[newFinalId]);
  combineShader.setUniform("front", GLSLValue(1));
  Shader::activeTexture(1);
  glBindTexture(GL_TEXTURE_2D, colorTexId[fullImgId]);

  combineShader.setUniform("front_depth", GLSLValue(2));
  Shader::activeTexture(2);
  glBindTexture(GL_TEXTURE_2D, depthTexId[fullImgId]);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_COMPARE_MODE, GL_NONE);

  combineShader.setUniform("back", GLSLValue(3));
  Shader::activeTexture(3);
  glBindTexture(GL_TEXTURE_2D, colorTexId[curPeelId]);

  combineShader.setUniform("back_depth", GLSLValue(4));
  Shader::activeTexture(4);
  glBindTexture(GL_TEXTURE_2D, depthTexId[curPeelId]);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_COMPARE_MODE, GL_NONE);

  combineShader.setUniform("volume1", GLSLValue(5));
  Shader::activeTexture(5);
  glBindTexture(GL_TEXTURE_2D, colorTexId[volume1]);

  combineShader.setUniform("volume2", GLSLValue(6));
  Shader::activeTexture(6);
  glBindTexture(GL_TEXTURE_2D, colorTexId[volume2]);

  Shader::activeTexture(Shader::AT_NONE);
  glBindTexture(GL_TEXTURE_2D, 0);

  combineShader.useShaders();
  combineShader.setupUniforms();

  drawColorTexture(-1);

  combineShader.stopUsingShaders();
  fullImgId = newFinalId;
}

void MorphoViewer::fastDraw()
{
  fast_draw = true;
  draw();
  fast_draw = false;
}

void MorphoViewer::draw() {
  draw(this);
}

void MorphoViewer::draw(QPaintDevice* device)
{
  if(quitting or !device)
    return;

  REPORT_GL_ERROR("draw: start");
  current_device = device;
  if(ImgData::MeshLineWidth == 0.0f) {
    GLfloat range[2];
    glGetFloatv(GL_LINE_WIDTH_RANGE, range);
    ImgData::MeshLineWidth = range[0];
  }
  updateFBOTex(drawWidth, drawHeight, fast_draw);

  uint oldSlices = ImgData::Slices;
  if(fast_draw)
    ImgData::Slices /= sampling;

  stack1->reloadLabelTex();
  stack2->reloadLabelTex();
  REPORT_GL_ERROR("draw: loading textures");

  // Set the clipping region
  // Draw the stack
  if(!raycastingShader1.initialized() or !raycastingShader2.initialized() or !volumeSurfShader1.initialized()
     or !volumeSurfShader2.initialized() or !textureSurfShader.initialized()) {
    Information::setStatus("Error, shaders not correctly initialized. Cannot render objects");
    return;
  }

  raycastingShader1.setUniform("tex", GLSLValue(Shader::AT_TEX3D));
  raycastingShader1.setUniform("colormap", GLSLValue(Shader::AT_CMAP_TEX));
  raycastingShader1.setUniform("second_tex", GLSLValue(Shader::AT_SECOND_TEX3D));
  raycastingShader1.setUniform("second_colormap", GLSLValue(Shader::AT_SECOND_CMAP_TEX));
  raycastingShader1.setUniform("front", GLSLValue(Shader::AT_FRONT_TEX));
  raycastingShader1.setUniform("front_color", GLSLValue(Shader::AT_FRONT_COLOR_TEX));
  raycastingShader1.setUniform("occlusion", GLSLValue(Shader::AT_OCCLUSION_TEX));
  raycastingShader1.setUniform("back", GLSLValue(Shader::AT_BACK_TEX));
  raycastingShader1.setUniform("solid", GLSLValue(Shader::AT_DEPTH_TEX));
  raycastingShader1.setUniform("labelcolormap", GLSLValue(Shader::AT_LABEL_TEX));
  raycastingShader1.setUniform("nb_colors", GLSLValue(int(ImgData::LabelColors.size())));

  raycastingShader2.setUniform("tex", GLSLValue(Shader::AT_TEX3D));
  raycastingShader2.setUniform("colormap", GLSLValue(Shader::AT_CMAP_TEX));
  raycastingShader2.setUniform("second_tex", GLSLValue(Shader::AT_SECOND_TEX3D));
  raycastingShader2.setUniform("second_colormap", GLSLValue(Shader::AT_SECOND_CMAP_TEX));
  raycastingShader2.setUniform("front", GLSLValue(Shader::AT_FRONT_TEX));
  raycastingShader2.setUniform("front_color", GLSLValue(Shader::AT_FRONT_COLOR_TEX));
  raycastingShader2.setUniform("occlusion", GLSLValue(Shader::AT_OCCLUSION_TEX));
  raycastingShader2.setUniform("back", GLSLValue(Shader::AT_BACK_TEX));
  raycastingShader2.setUniform("solid", GLSLValue(Shader::AT_DEPTH_TEX));
  raycastingShader2.setUniform("labelcolormap", GLSLValue(Shader::AT_LABEL_TEX));
  raycastingShader2.setUniform("nb_colors", GLSLValue(int(ImgData::LabelColors.size())));

  volumeSurfShader1.setUniform("backDepth", GLSLValue(Shader::AT_DEPTH_TEX));
  volumeSurfShader1.setUniform("frontDepth", GLSLValue(Shader::AT_FRONT_TEX));
  volumeSurfShader1.setUniform("testFront", GLSLValue(false));
  volumeSurfShader1.setUniform("tex", GLSLValue(Shader::AT_TEX3D));
  volumeSurfShader1.setUniform("colormap", GLSLValue(Shader::AT_CMAP_TEX));
  volumeSurfShader1.setUniform("second_tex", GLSLValue(Shader::AT_SECOND_TEX3D));
  volumeSurfShader1.setUniform("second_colormap", GLSLValue(Shader::AT_SECOND_CMAP_TEX));
  volumeSurfShader1.setUniform("labelcolormap", GLSLValue(Shader::AT_LABEL_TEX));
  volumeSurfShader1.setUniform("nb_colors", GLSLValue(int(ImgData::LabelColors.size())));

  volumeSurfShader2.setUniform("backDepth", GLSLValue(Shader::AT_DEPTH_TEX));
  volumeSurfShader2.setUniform("frontDepth", GLSLValue(Shader::AT_FRONT_TEX));
  volumeSurfShader2.setUniform("testFront", GLSLValue(false));
  volumeSurfShader2.setUniform("tex", GLSLValue(Shader::AT_TEX3D));
  volumeSurfShader2.setUniform("colormap", GLSLValue(Shader::AT_CMAP_TEX));
  volumeSurfShader2.setUniform("second_tex", GLSLValue(Shader::AT_SECOND_TEX3D));
  volumeSurfShader2.setUniform("second_colormap", GLSLValue(Shader::AT_SECOND_CMAP_TEX));
  volumeSurfShader2.setUniform("labelcolormap", GLSLValue(Shader::AT_LABEL_TEX));
  volumeSurfShader2.setUniform("nb_colors", GLSLValue(int(ImgData::LabelColors.size())));

  indexSurfShader.setUniform("backDepth", GLSLValue(Shader::AT_DEPTH_TEX));
  indexSurfShader.setUniform("frontDepth", GLSLValue(Shader::AT_FRONT_TEX));
  indexSurfShader.setUniform("testFront", GLSLValue(false));
  indexSurfShader.setUniform("surfcolormap", GLSLValue(Shader::AT_SURF_TEX));
  indexSurfShader.setUniform("labelcolormap", GLSLValue(Shader::AT_LABEL_TEX));
  indexSurfShader.setUniform("heatcolormap", GLSLValue(Shader::AT_HEAT_TEX));
  indexSurfShader.setUniform("nb_colors", GLSLValue(int(ImgData::LabelColors.size())));
  indexSurfShader.setUniform("labels", GLSLValue(false));
  indexSurfShader.setUniform("heatmap", GLSLValue(false));

  textureSurfShader.setUniform("backDepth", GLSLValue(Shader::AT_DEPTH_TEX));
  textureSurfShader.setUniform("frontDepth", GLSLValue(Shader::AT_FRONT_TEX));
  textureSurfShader.setUniform("testFront", GLSLValue(false));

  REPORT_GL_ERROR("draw: setting up shaders");

  // Save projection and modelview matrices
  Color3f clearColor = Colors::getColor(Colors::BackgroundColor);
  glClearColor(clearColor.r(), clearColor.g(), clearColor.b(), 1.0);
  REPORT_GL_ERROR("draw: Set ClearColor");

  setupFramebuffer(depthTexId[FI_BACKGROUND], colorTexId[FI_BACKGROUND]);
  REPORT_GL_ERROR("draw: Setting up FrameBuffer");

  textureSurfShader.setUniform("peeling", GLSLValue(false));
  indexSurfShader.setUniform("peeling", GLSLValue(false));
  volumeSurfShader1.setUniform("peeling", GLSLValue(false));
  volumeSurfShader2.setUniform("peeling", GLSLValue(false));
  REPORT_GL_ERROR("draw: Turn off surface shader peeling");

  // ### Begin draw opaque

  glDisable(GL_BLEND);

  clip1.drawGrid(SceneRadius);
  clip2.drawGrid(SceneRadius);
  clip3.drawGrid(SceneRadius);
  REPORT_GL_ERROR("draw: clipping plane grid");

  if(gridIsDrawn()) {
    glLineWidth(1.0);
    drawGrid(camera()->sceneRadius());
    REPORT_GL_ERROR("draw: grid");
  }
	if(axisIsDrawn()) {
    glLineWidth(2.0);
    drawAxis(camera()->sceneRadius());
    REPORT_GL_ERROR("draw: axis");
  }

  stack1->drawBBox();
  stack2->drawBBox();
  REPORT_GL_ERROR("draw: BBox");

  clip1.drawClip();
  clip2.drawClip();
  clip3.drawClip();
  REPORT_GL_ERROR("draw: clipping planes");

  stack1->drawMesh();
  stack2->drawMesh();
  REPORT_GL_ERROR("draw: drawing mesh");

  glEnable(GL_LIGHTING);

  if(stack1->showOpaqueSurface()) {
    setLighting(stack1);
    stack1->drawSurf(false, 
        &textureSurfShader, &indexSurfShader, &colorSurfShader, &volumeSurfShader1);
  }
  if(stack2->showOpaqueSurface()) {
    setLighting(stack2);
    stack2->drawSurf(false, 
        &textureSurfShader, &indexSurfShader, &colorSurfShader, &volumeSurfShader2);
  }
  REPORT_GL_ERROR("draw: drawing surface");

  cutSurf->drawCutSurfGrid(*stack1);
  cutSurf->drawCutSurfGrid(*stack2);

  if(cutSurf->showOpaqueSurface(*stack1)) {
    setLighting(stack1);
    cutSurf->drawCutSurf(*stack1, false, &volumeSurfShader1);
  }

  if(cutSurf->showOpaqueSurface(*stack2)) {
    setLighting(stack2);
    cutSurf->drawCutSurf(*stack2, false, &volumeSurfShader2);
  }
  REPORT_GL_ERROR("draw: drawing cutting surface");

  stack1->drawAxis();
  stack2->drawAxis();
  stack1->drawVertexVertexLine(stack2);
  stack2->drawVertexVertexLine(stack1);
  REPORT_GL_ERROR("draw: axis");

  // ### End draw opaque

  glClearColor(0, 0, 0, 0);

  Shader::activeTexture(Shader::AT_DEPTH_TEX);
  glBindTexture(GL_TEXTURE_2D, depthTexId[FI_BACKGROUND]);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_COMPARE_MODE, GL_NONE);
  Shader::activeTexture(Shader::AT_NONE);

  glEnable(GL_DEPTH_TEST);
  REPORT_GL_ERROR("draw: setting 2D texture");

  int fullImgId = FI_FULL_IMG1;

  // First, clean the first peel and set it "full front"
  glClearDepth(0.0);
  glClearColor(0, 0, 0, 0);
  resetupFramebuffer(depthTexId[fullImgId], colorTexId[fullImgId]);

  glClearDepth(1.0);

  textureSurfShader.setUniform("peeling", GLSLValue(true));
  indexSurfShader.setUniform("peeling", GLSLValue(true));
  volumeSurfShader1.setUniform("peeling", GLSLValue(true));
  volumeSurfShader2.setUniform("peeling", GLSLValue(true));

  // Clear occlusion texture
  glClearColor(0, 1, 0, 1);
  resetupFramebuffer(0, colorTexId[FI_OCCLUSION]);

  glClearColor(0, 0, 0, 0);
  REPORT_GL_ERROR("draw: setting 3D draw shaders");

  // Create query
  GLuint occ_query, nb_pixels;
  glfuncs->glGenQueries(1, &occ_query);
  for(int current_pass = 0; current_pass < MaxNbPeels; ++current_pass) {
    // Setup the last full image as front texture

    // Set the framebuffer to write the current peel

    resetupFramebuffer(depthTexId[FI_CUR_PEEL], colorTexId[FI_CUR_PEEL]);

    if(current_pass == 0 and !fast_draw) {
      // Now, copy the depth buffer, only on the first pass, to get the correct
      // depth of the surface
      setupCopyFB(depthTexId[FI_BACKGROUND], colorTexId[FI_BACKGROUND]);

      glfuncs->glBindFramebuffer(GL_READ_FRAMEBUFFER_EXT, fboCopyId);
      glfuncs->glBindFramebuffer(GL_DRAW_FRAMEBUFFER_EXT, fboId);

      glfuncs->glBlitFramebuffer(0, 0, texWidth, texHeight, 0, 0, texWidth, texHeight, GL_DEPTH_BUFFER_BIT, GL_NEAREST);

      glfuncs->glBindFramebuffer(GL_FRAMEBUFFER_EXT, fboId);
    }

    Shader::activeTexture(Shader::AT_DEPTH_TEX);
    glBindTexture(GL_TEXTURE_2D, depthTexId[FI_BACKGROUND]);
    Shader::activeTexture(Shader::AT_FRONT_TEX);
    glBindTexture(GL_TEXTURE_2D, depthTexId[fullImgId]);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_COMPARE_MODE, GL_NONE);
    Shader::activeTexture(Shader::AT_NONE);

    // ### Begin draw transparent

    glEnable(GL_LIGHTING);
    if(stack1->showTransparentSurface()) {
      setLighting(stack1);
      stack1->drawSurf(false, &textureSurfShader, &indexSurfShader, &volumeSurfShader1);
    }
    if(stack2->showTransparentSurface()) {
      setLighting(stack2);
      stack2->drawSurf(false, &textureSurfShader, &indexSurfShader, &volumeSurfShader2);
    }
    // Next, the cutting surface
    if(cutSurf->showTransparentSurface(*stack1)) {
      setLighting(stack1);
      cutSurf->drawCutSurf(*stack1, false, &volumeSurfShader1);
    }
    if(cutSurf->showTransparentSurface(*stack1)) {
      setLighting(stack2);
      cutSurf->drawCutSurf(*stack2, false, &volumeSurfShader2);
    }

    // ### End draw transparent

    // And the volume

    // Setup depth buffers
    Shader::activeTexture(Shader::AT_BACK_TEX);
    glBindTexture(GL_TEXTURE_2D, depthTexId[FI_CUR_PEEL]);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_COMPARE_MODE, GL_NONE);
    Shader::activeTexture(Shader::AT_FRONT_TEX);
    glBindTexture(GL_TEXTURE_2D, depthTexId[fullImgId]);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_COMPARE_MODE, GL_NONE);
    Shader::activeTexture(Shader::AT_FRONT_COLOR_TEX);
    glBindTexture(GL_TEXTURE_2D, colorTexId[fullImgId]);
    Shader::activeTexture(Shader::AT_OCCLUSION_TEX);
    glBindTexture(GL_TEXTURE_2D, colorTexId[FI_OCCLUSION]);
    Shader::activeTexture(Shader::AT_NONE);
    glBindTexture(GL_TEXTURE_2D, 0);

    resetupFramebuffer(0, colorTexId[FI_VOLUME1]);
    stack1->drawStack(&raycastingShader1);

    resetupFramebuffer(0, colorTexId[FI_VOLUME2]);
    stack2->drawStack(&raycastingShader2);

    combinePeels(fullImgId, FI_CUR_PEEL, FI_VOLUME1, FI_VOLUME2);

    // Check if there is anything else to be done

    resetupFramebuffer(0, colorTexId[FI_OCCLUSION]);

    Shader::activeTexture(1);
    glBindTexture(GL_TEXTURE_2D, depthTexId[fullImgId]);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_COMPARE_MODE, GL_NONE);
    Shader::activeTexture(2);
    glBindTexture(GL_TEXTURE_2D, colorTexId[fullImgId]);
    Shader::activeTexture(3);
    glBindTexture(GL_TEXTURE_2D, depthTexId[FI_BACKGROUND]);
    Shader::activeTexture(Shader::AT_NONE);
    glBindTexture(GL_TEXTURE_2D, 0);

    occlusionShader.setUniform("depth", GLSLValue(1));
    occlusionShader.setUniform("color", GLSLValue(2));
    occlusionShader.setUniform("background", GLSLValue(3));

    occlusionShader.useShaders();
    occlusionShader.setupUniforms();

    glfuncs->glBeginQuery(GL_SAMPLES_PASSED_ARB, occ_query);

    drawColorTexture(-1);

    glfuncs->glEndQuery(GL_SAMPLES_PASSED_ARB);

    glfuncs->glGetQueryObjectuiv(occ_query, GL_QUERY_RESULT_ARB, &nb_pixels);

    occlusionShader.stopUsingShaders();

    if(DEBUG and current_pass == MaxNbPeels - 1 and nb_pixels > 0) {
      int prevFullImg = fullImgId;
      fullImgId = (fullImgId == FI_FULL_IMG1) ? FI_FULL_IMG2 : FI_FULL_IMG1;
      resetupFramebuffer(0, colorTexId[fullImgId]);
      glDisable(GL_DEPTH_TEST);
      drawColorTexture(prevFullImg);
      glEnable(GL_BLEND);
      drawColorTexture(FI_OCCLUSION);
      glDisable(GL_BLEND);
      glEnable(GL_DEPTH_TEST);
      Information::setStatus(QString("Number of left pixels on pass # %1 = %2.").arg(current_pass).arg(nb_pixels));
    }

    // End check

    if(current_pass == 0 and !fast_draw) {
      // Now, copy the depth buffer!
      glEnable(GL_TEXTURE_2D);
      glBindTexture(GL_TEXTURE_2D, depthTexId[FI_CUR_PEEL]);
      depthTexture.resize(texWidth * texHeight);
      glGetTexImage(GL_TEXTURE_2D, 0, GL_DEPTH_COMPONENT, GL_FLOAT, &depthTexture[0]);
      glBindTexture(GL_TEXTURE_2D, 0);
      glDisable(GL_TEXTURE_2D);
    }
    if(show_slice > -1) {
      if(current_pass == show_slice) {
        int prevFullImg = fullImgId;
        fullImgId = (fullImgId == FI_FULL_IMG1) ? FI_FULL_IMG2 : FI_FULL_IMG1;

        glBindTexture(GL_TEXTURE_2D, depthTexId[FI_CUR_PEEL]);
        glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_COMPARE_MODE, GL_NONE);
        glBindTexture(GL_TEXTURE_2D, depthTexId[prevFullImg]);
        glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_COMPARE_MODE, GL_NONE);
        glBindTexture(GL_TEXTURE_2D, 0);

        // glEnable(GL_BLEND);
        // glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
        glDisable(GL_DEPTH_TEST);

        glClearColor(0, 0, 0, 1);
        resetupFramebuffer(0, colorTexId[fullImgId]);
        switch(slice_type) {
        case 0:
          drawColorTexture(FI_CUR_PEEL);
          break;
        case 1:
          drawColorTexture(FI_CUR_PEEL, true);
          break;
        case 2:
          drawColorTexture(FI_VOLUME1);
          break;
        case 3:
          drawColorTexture(FI_VOLUME2);
          break;
        case 4:
          drawColorTexture(prevFullImg);
          break;
        case 5:
          drawColorTexture(prevFullImg, true);
          break;
        case 6:
          drawColorTexture(FI_OCCLUSION);
          break;
        default:
          Information::out << "Unknown slice type: " << slice_type << endl;
        }
        break;
      }
    }
    REPORT_GL_ERROR("draw");
    if(nb_pixels == 0)     // TODO: Find out why some pixels are never empty!
    {
      if(DEBUG)
        Information::out << "# peels = " << current_pass + 1 << endl;
      if(show_slice > current_pass) {
        int prevFullImg = fullImgId;
        fullImgId = (fullImgId == FI_FULL_IMG1) ? FI_FULL_IMG2 : FI_FULL_IMG1;
        glDisable(GL_DEPTH_TEST);

        glBindTexture(GL_TEXTURE_2D, depthTexId[prevFullImg]);
        glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_COMPARE_MODE, GL_NONE);

        glClearColor(0, 0, 0, 1);
        resetupFramebuffer(0, colorTexId[fullImgId]);
        switch(slice_type) {
        case 0:
          drawColorTexture(FI_BACKGROUND);
          break;
        case 1:
          drawColorTexture(FI_BACKGROUND, true);
          break;
        case 5:
          drawColorTexture(prevFullImg, true);
          break;
        default:
          drawColorTexture(prevFullImg);
          break;
        }
        break;
      }
      break;
    } else {
      if(DEBUG)
        Information::out << "# pixels at step " << current_pass << " = " << nb_pixels << endl;
    }
  }
  REPORT_GL_ERROR("draw: drawing 3D volume");
  glfuncs->glDeleteQueries(1, &occ_query);

  int finalId = (fullImgId == FI_FULL_IMG1) ? FI_FULL_IMG2 : FI_FULL_IMG1;
  resetupFramebuffer(0, colorTexId[finalId]);
  REPORT_GL_ERROR("draw: reseting frame buffer");

  Shader::activeTexture(Shader::AT_DEPTH_TEX);
  glBindTexture(GL_TEXTURE_2D, 0);
  Shader::activeTexture(Shader::AT_NONE);
  REPORT_GL_ERROR("draw: reseting 2D texture");

  // Now, render on the final texture, before the screen

  clip1.disable();
  clip2.disable();
  clip3.disable();

  // drawColorTexture(0);

  {
    resetupFramebuffer(depthTexId[finalId], colorTexId[finalId]);
    finalCombineShader.setUniform("front", GLSLValue(1));
    Shader::activeTexture(1);
    glBindTexture(GL_TEXTURE_2D, colorTexId[fullImgId]);

    finalCombineShader.setUniform("back", GLSLValue(3));
    Shader::activeTexture(3);
    glBindTexture(GL_TEXTURE_2D, colorTexId[FI_BACKGROUND]);

    Shader::activeTexture(Shader::AT_NONE);
    glBindTexture(GL_TEXTURE_2D, 0);

    finalCombineShader.useShaders();
    finalCombineShader.setupUniforms();

    drawColorTexture(-1);

    finalCombineShader.stopUsingShaders();
  }
  REPORT_GL_ERROR("draw: rendering to final texture");

  /*
     glEnable(GL_BLEND);
     glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
     glDisable(GL_DEPTH_TEST);

     if(show_slice == -1)
     drawColorTexture(FI_BACKGROUND);

     if(fullImgId >= 0)
     {
     drawColorTexture(fullImgId);
     }
   */

  resetFramebuffer();
  REPORT_GL_ERROR("draw");

  glDisable(GL_DEPTH_TEST);

  postProcessShader.activeTexture(Shader::AT_FINAL_VOLUME_TEX);
  glBindTexture(GL_TEXTURE_2D, colorTexId[finalId]);
  postProcessShader.activeTexture(Shader::AT_NONE);
  postProcessShader.setUniform("texId", GLSLValue(Shader::AT_FINAL_VOLUME_TEX));
  postProcessShader.setUniform("texSize", GLSLValue(Point2i(texWidth, texHeight)));
  postProcessShader.setUniform("brightness", GLSLValue(GlobalBrightness));
  postProcessShader.setUniform("contrast", GLSLValue(GlobalContrast));
  /*
   * if(not fast_draw and UnsharpStrength != 0)
   *{
   *  post_processShader.setUniform("unsharp", GLSLValue(true));
   *  post_processShader.setUniform("amount", GLSLValue(UnsharpStrength));
   *  post_processShader.setUniform("kernel", GLSLValue(unsharp_kernel, 9));
   *}
   * else
   */
  postProcessShader.setUniform("unsharp", GLSLValue(false));
  postProcessShader.useShaders();
  postProcessShader.setupUniforms();
  REPORT_GL_ERROR("draw: post processing shaders");

  drawColorTexture(-1);
  REPORT_GL_ERROR("draw: final draw");

  postProcessShader.stopUsingShaders();


  // ### Begin draw overlay

  // Draw cell correspondence
  drawCellMap();

  drawSelectRect();
  drawSelectLasso();
  drawVoxelCursor();
  REPORT_GL_ERROR("draw: selection tools");

  // Draw the legend and scale bar
  drawColorBar();
  drawScaleBar();
  REPORT_GL_ERROR("draw: color bar");

  // ### End draw overlay

  glEnable(GL_DEPTH_TEST);

  if(fast_draw)
    ImgData::Slices = oldSlices;

  REPORT_GL_ERROR("draw: final");
}

void MorphoViewer::clipEnable()
{
  // Turn on clipping planes
  clip1.drawClip();
  clip2.drawClip();
  clip3.drawClip();
}

void MorphoViewer::clipDisable()
{
  // Turn off clipping planes
  clip1.disable();
  clip2.disable();
  clip3.disable();
}

void MorphoViewer::drawScaleBar()
{
  try {
    ImgData::scaleBar.setScale(1e-6); // We work in microns
    ImgData::scaleBar.draw(this, current_device);
  }
  catch(const QString& s) {
    Information::out << "MorphoViewer::drawScaleBar() - " << s << endl;
  }
  catch(...) {
    Information::out << "MorphoViewer::drawScaleBar() - Unknown exception" << endl;
  }
}

void MorphoViewer::drawColorBar()
{
  Mesh &mesh1 = *stack1->mesh;
  bool stk1 = mesh1.showSignal() or mesh1.showHeat() or mesh1.showWallHeat() or mesh1.showTriangleValue();

  Mesh &mesh2 = *stack2->mesh;
  bool stk2 = mesh2.showSignal() or mesh2.showHeat() or mesh2.showWallHeat() or mesh2.showTriangleValue();

  // If no surfaces, or both surfaces and selection different return
  // Note that you can still different heat maps on the two stacks, so this is not really correct
  // Best would be to display both Colorbars.
  if(!stk1 and !stk2) {
    return;
  } else if(stk1 and stk2 and 
           (mesh1.showSignal() != mesh2.showSignal() or mesh1.showHeat() != mesh2.showHeat() 
             or mesh1.showWallHeat() != mesh2.showWallHeat() or mesh1.showTriangleValue() != mesh2.showTriangleValue()))
    return;

  // Get the stack and view
  ImgData* stk = (stk1 ? stack1 : stack2);

  if(stk->mesh->showHeat() or stk->mesh->showWallHeat()) {
    stk->setupHeatColorMap();
    ImgData::colorBar.vmin = stk->mesh->heatMapBounds()[0];
    ImgData::colorBar.vmax = stk->mesh->heatMapBounds()[1];
    ImgData::colorBar.label = stk->mesh->heatMapUnit();
    ImgData::colorBar.draw(stk->heatTexId, current_device);
  } else if(stk->mesh->showSignal()) {
    stk->setupSurfColorMap();
    ImgData::colorBar.vmin = stk->mesh->signalBounds()[0];
    ImgData::colorBar.vmax = stk->mesh->signalBounds()[1];
    ImgData::colorBar.label = stk->mesh->signalUnit();
    ImgData::colorBar.draw(stk->surfTexId, current_device);
  } else if(stk->mesh->showTriangleValue()) {
    stk->setupTexFromColorMap(stk->colorMapTexId, stk->mesh->triangleValueColors());
    ImgData::colorBar.vmin = stk->mesh->triangleValueBounds()[0];
    ImgData::colorBar.vmax = stk->mesh->triangleValueBounds()[1];
    ImgData::colorBar.label = stk->mesh->triangleValueUnit();
    ImgData::colorBar.draw(stk->colorMapTexId, current_device);
  }
}

void MorphoViewer::drawCellMap()
{
  if(!mainWindow()->runningProcess()) {
    bool stack1Draw = (stack1->mesh->showMeshCellMap() 
                       and stack1->pntsVA.size() > 0 and !stack1->mesh->labelCenter().empty());
    bool stack2Draw = (stack2->mesh->showMeshCellMap() 
                       and stack2->pntsVA.size() > 0 and !stack2->mesh->labelCenter().empty());

    if(!stack1Draw and !stack2Draw)
      return;

    glDisable(GL_TEXTURE_1D);
    glDisable(GL_TEXTURE_2D);
    glDisable(GL_TEXTURE_3D);
    glDisable(GL_BLEND);
    glDisable(GL_AUTO_NORMAL);
    glDisable(GL_LIGHTING);

    // glPolygonMode(GL_FRONT, GL_LINE);
    // glLineWidth(ImgData::MeshLineWidth);
    // Color3f meshColor = Colors::getColor(stack1->MeshColor);
    // Color3f meshColor = ImgData::palette->getColor(stack1->MeshColor);
    // meshColor /= 255.0;
    // glColor3fv(meshColor.data());

    // if(stack1Draw and stack2Draw) {
    //  glBegin(GL_LINES);
    //  forall(const IntPoint3fPair &p, stack1->LabelCenter)
    //    if(p.first > 0 and stack1->LabelCenter.find(p.first) != stack1->LabelCenter.end()
    //                 and stack2->LabelCenter.find(p.first) != stack2->LabelCenter.end()) {
    //      Point3f p1 = stack1->LabelCenter[p.first];
    //      Point3f p2 = stack2->LabelCenter[p.first];
    //      p1 = Point3f(stack1->getFrame()->inverseCoordinatesOf(Vec(p1)));
    //      p2 = Point3f(stack2->getFrame()->inverseCoordinatesOf(Vec(p2)));
    //      glVertex3fv(p1.c_data());
    //      glVertex3fv(p2.c_data());
    //    }
    //  glEnd();
    //}
    glDisable(GL_DEPTH_TEST);
    if(stack1Draw) {
      glColor3fv(Colors::getColor(Colors::Mesh1CellsColor).c_data());
      forall(const IntPoint3fPair& p, stack1->mesh->labelCenter()) {
        Point3f pos = Point3f(camera()->projectedCoordinatesOf(Vec(p.second), &stack1->getFrame()));
        // Information::out << "Pos:" << pos << endl;
        drawText((int)pos[0], (int)pos[1], QString::number(p.first));
      }
    }
    if(stack2Draw) {
      glColor3fv(Colors::getColor(Colors::Mesh2CellsColor).c_data());
      forall(const IntPoint3fPair& p, stack2->mesh->labelCenter()) {
        Point3f pos = Point3f(camera()->projectedCoordinatesOf(Vec(p.second), &stack2->getFrame()));
        drawText((int)pos[0], (int)pos[1], QString::number(p.first));
      }
    }
  }
}

void MorphoViewer::drawVoxelCursor()
{
  if(!voxelEditCursor)
    return;

  {
    startScreenCoordinatesSystem();
    glDisable(GL_LIGHTING);
    glDisable(GL_TEXTURE_1D);
    glDisable(GL_TEXTURE_2D);
    glDisable(GL_TEXTURE_3D);

    Color3f pixelColor = Colors::getColor(Colors::VoxelEditColor);

    float voxelEditRadius = float(ImgData::VoxelEditRadius);

    if(DEBUG) {
      Information::out << "Drawing voxelEditCursor around mouse position " 
         << mousePos.x() << "x" << mousePos.y() << endl;
      Information::out << "   Color: " << pixelColor << " - Radius: " << voxelEditRadius << endl;
    }

    glLineWidth(1.0);
    Point2f c(mousePos.x(), mousePos.y());
    glBegin(GL_LINE_LOOP);
    glColor4fv(pixelColor.c_data());
    glNormal3f(0, 0, 1);
    for(int i = 0; i < 24; ++i) {
      double a = 2 * M_PI * float(i) / 24;
      Point2f p = c + voxelEditRadius * Point2f(cos(a), sin(a));
      glVertex2fv(p.c_data());
    }
    glEnd();
    stopScreenCoordinatesSystem();
  }

  /*
     {
     QPainter paint(this);
     paint.setRenderHint(QPainter::HighQualityAntialiasing);
     QPen pen(Colors::getColor(Colors::VoxelEditColor));
     pen.setWidth(0);
     paint.setPen(pen);
     paint.setBrush(Qt::NoBrush);
     paint.drawEllipse(mousePos, ImgData::VoxelEditRadius, ImgData::VoxelEditRadius);
     }
   */

  if(DrawClipBox) {
    glDisable(GL_LIGHTING);
    glDisable(GL_TEXTURE_1D);
    glDisable(GL_TEXTURE_2D);
    glDisable(GL_TEXTURE_3D);

    ImgData* stk = findSelectStack();
    if(stk) {
      Point3f bBox[2] = { stk->imageToWorld(stk->bBoxTex[0]), stk->imageToWorld(stk->bBoxTex[1]) };
      if(bBox[0].x() < bBox[1].x() and bBox[0].y() < bBox[1].y() and bBox[0].z() < bBox[1].z()) {
        glMatrixMode(GL_MODELVIEW);
        glPushMatrix();
        glMultMatrixd(stk->getFrame().worldMatrix());
        glBegin(GL_LINES);
        for(int i = 0; i < 2; i++)
          for(int j = 0; j < 2; j++) {
            for(int k = 0; k < 2; k++)
              glVertex3f(bBox[k].x(), bBox[i].y(), bBox[j].z());
            for(int k = 0; k < 2; k++)
              glVertex3f(bBox[i].x(), bBox[k].y(), bBox[j].z());
            for(int k = 0; k < 2; k++)
              glVertex3f(bBox[i].x(), bBox[j].y(), bBox[k].z());
          }
        glEnd();

        glPopMatrix();
      }
    }
  }
}

void MorphoViewer::drawSelectRect()
{
  if(guiActionOn and guiAction == MESH_SEL_POINTS and selectRect.topLeft() == selectRect.bottomRight())
    return;

  startScreenCoordinatesSystem();
  glDisable(GL_LIGHTING);
  glDisable(GL_TEXTURE_1D);
  glDisable(GL_TEXTURE_2D);
  glDisable(GL_TEXTURE_3D);
  glEnable(GL_BLEND);
  glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
  glDisable(GL_DEPTH_TEST);

  glColor4f(0.2, 0.2, 0.2f, 0.5f);
  glPolygonMode(GL_FRONT, GL_FILL);
  glBegin(GL_QUADS);
  glVertex2i(selectRect.left(), selectRect.top());
  glVertex2i(selectRect.right(), selectRect.top());
  glVertex2i(selectRect.right(), selectRect.bottom());
  glVertex2i(selectRect.left(), selectRect.bottom());
  glEnd();

  glLineWidth(2.0);
  glColor4f(0.5f, 0.5f, 0.5f, 0.5f);
  glBegin(GL_LINE_LOOP);
  glVertex2i(selectRect.left(), selectRect.top());
  glVertex2i(selectRect.right(), selectRect.top());
  glVertex2i(selectRect.right(), selectRect.bottom());
  glVertex2i(selectRect.left(), selectRect.bottom());
  glEnd();

  stopScreenCoordinatesSystem();
}

void MorphoViewer::drawSelectLasso()
{
  if(guiActionOn and guiAction == MESH_SEL_POINTS 
                            and selectLasso.contour.size()==0)
    return;

  startScreenCoordinatesSystem();
  glDisable(GL_LIGHTING);
  glDisable(GL_TEXTURE_1D);
  glDisable(GL_TEXTURE_2D);
  glDisable(GL_TEXTURE_3D);
  glEnable(GL_BLEND);
  glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
  glDisable(GL_DEPTH_TEST);

  glLineWidth(3.0);
  glColor4f(0.7f, 0.7f, 0.7f, 0.7f);
  glBegin(GL_LINE_LOOP);
  for(uint i=0;i<selectLasso.contour.size();i++)
  glVertex2i(selectLasso.contour[i].x(),selectLasso.contour[i].y());
  glEnd();

  stopScreenCoordinatesSystem();
}

void MorphoViewer::setLighting(ImgData *stack)
{
  glEnable(GL_COLOR_MATERIAL);
  glColorMaterial(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE);
  glMatrixMode(GL_MODELVIEW);
  glPushMatrix();
  glLoadIdentity();

  // Set lighting, 4 lights
  // GLfloat lightc[4] = {stack->mesh->brightness()/8.0f, stack->mesh->brightness()/8.0f,
  // stack->mesh->brightness()/8.0f, 1.0f };
  float br = stack->mesh->brightness();
  float lc = br / 8.0f;
  float la = br / 8.0f;
  float ls = GlobalSpecular * br;
  Colorf lightc(lc, lc, lc, 1.0f);
  Colorf lights(ls, ls, ls, 1.0f);
  Colorf lighta(la, la, la, 1.0f);

  for(int i = 0; i < 4; ++i) {
    glLightfv(GL_LIGHT0 + i, GL_AMBIENT, lighta.c_data());
    glLightfv(GL_LIGHT0 + i, GL_DIFFUSE, lightc.c_data());
    glLightfv(GL_LIGHT0 + i, GL_SPECULAR, lights.c_data());
  }

  glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, GlobalShininess);
  glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, lights.c_data());

  Point4f lightp0(1.0f, -1.0f, 1.0f, 0.0f);
  Point4f lightp1(-1.0f, -1.0f, 1.0f, 0.0f);
  Point4f lightp2(1.0f, 1.0f, 1.0f, 0.0f);
  Point4f lightp3(-1.0f, 1.0f, 1.0f, 0.0f);

  glLightfv(GL_LIGHT0, GL_POSITION, lightp0.c_data());
  glLightfv(GL_LIGHT1, GL_POSITION, lightp1.c_data());
  glLightfv(GL_LIGHT2, GL_POSITION, lightp2.c_data());
  glLightfv(GL_LIGHT3, GL_POSITION, lightp3.c_data());

  glEnable(GL_LIGHT0);
  glEnable(GL_LIGHT1);
  glEnable(GL_LIGHT2);
  glEnable(GL_LIGHT3);

  glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, GL_TRUE);
  Colorf black(0, 0, 0, 1);
  glLightModelfv(GL_LIGHT_MODEL_AMBIENT, black.c_data());

  glPopMatrix();
}

// Read from parameter file
void MorphoViewer::readParms(Parms& parms, QString section)
{
  parms(section, "DrawCellMap", DrawCellMap, false);
  parms(section, "DrawClipBox", DrawClipBox, false);
  parms(section, "Spinning", Spinning, 10000.0f);
  camera()->frame()->setSpinningSensitivity(Spinning);
  parms(section, "FlySpeed", FlySpeed, .005f);
  parms(section, "Brightness", GlobalBrightness, .0f);
  parms(section, "Contrast", GlobalContrast, 1.0f);
  parms(section, "Shininess", GlobalShininess, 32.0f);
  parms(section, "Specular", GlobalSpecular, 0.2f);
  parms(section, "UnsharpStrength", UnsharpStrength, 1.0f);
  parms(section, "MaxNbPeels", MaxNbPeels, 10);
  if(MaxNbPeels < 1)
    MaxNbPeels = 1;
  CameraFrame = Matrix4d::identity();;
  CameraFrame[3][2] = -1000;
  parms(section, "CameraFrame", CameraFrame, CameraFrame);
  camera()->setFromModelViewMatrix(CameraFrame.data());
  parms(section, "SceneRadius", SceneRadius, 500.0f);
  camera()->setSceneRadius(SceneRadius);

  float zoom;
  parms(section, "CameraZoom", zoom, 1.0f);
  _camera->setZoom(zoom);

  //camera()->showEntireScene();

  clip1.readParms(parms, "Clip1");
  clip2.readParms(parms, "Clip2");
  clip3.readParms(parms, "Clip3");
}

// Write to parameter file
void MorphoViewer::writeParms(QTextStream& pout, QString section)
{
  pout << endl;
  pout << "[" << section << "]" << endl;
  pout << "DrawCellMap: " << (DrawCellMap ? "true" : "false") << endl;
  pout << "DrawClipBox: " << (DrawClipBox ? "true" : "false") << endl;
  pout << "Spinning: " << Spinning << endl;
  pout << "FlySpeed: " << FlySpeed << endl;
  pout << "Brightness: " << GlobalBrightness << endl;
  pout << "Contrast: " << GlobalContrast << endl;
  pout << "Shininess: " << GlobalShininess << endl;
  pout << "Specular: " << GlobalSpecular << endl;
  pout << "UnsharpStrength: " << UnsharpStrength << endl;
  pout << "MaxNbPeels: " << MaxNbPeels << endl;
  camera()->getModelViewMatrix(CameraFrame.data());
  pout << "CameraFrame: " << CameraFrame << endl;
  pout << "SceneRadius: " << SceneRadius << endl;
  pout << "CameraZoom: " << _camera->zoom() << endl;
  pout << endl;

  clip1.writeParms(pout, "Clip1");
  clip2.writeParms(pout, "Clip2");
  clip3.writeParms(pout, "Clip3");
}

// Set the current label and color
void MorphoViewer::setLabel(int label)
{
  uint idx = label % ImgData::LabelColors.size();
  if(idx >= ImgData::LabelColors.size())
    return;

  QPixmap pix(24, 24);
  if(label == 0)
    pix.fill(QColor(0, 0, 0, 0));
  else {
    Point4f col = Point4f(ImgData::LabelColors[idx]);
    pix.fill(QColor(col.x() * 255, col.y() * 255, col.z() * 255));
  }
  QIcon icon = QIcon(pix);
  selectedLabel = label;
  emit selectLabelChanged(label);
  emit setLabelColor(icon);
}

MorphoGraphX* MorphoViewer::mainWindow()
{
  QWidget* p = this;
  MorphoGraphX* mgx = 0;
  do {
    p = p->parentWidget();
    mgx = dynamic_cast<MorphoGraphX*>(p);
  } while(p and !mgx);
  return mgx;
}

// Find which stack to use for pixel editing
ImgData* MorphoViewer::findSelectStack()
{
  if(processRunning())
    return 0;
  MorphoGraphX* mgx = mainWindow();
  if(!mgx) {
    Information::out << "Could not find main window" << endl;
    return 0;
  }
  if(mgx->activeStack() == 0 and stack1->stack->work()->isVisible() and stack1->valid())
    return stack1;
  if(mgx->activeStack() == 1 and stack2->stack->work()->isVisible() and stack2->valid())
    return stack2;
  return (0);
}

// Find which stack to use for label selection
ImgData* MorphoViewer::findSelectSurf()
{
  if(processRunning())
    return 0;
  MorphoGraphX* mgx = mainWindow();
  if(!mgx) {
    Information::out << "Could not find main window" << endl;
    return 0;
  }
  if(mgx->activeMesh() == 0 and stack1->mesh->showSurface())   // and stack1->SurfLabels)
    return stack1;
  if(mgx->activeMesh() == 1 and stack2->mesh->showSurface())   // and stack2->SurfLabels)
    return stack2;
  return (0);
}

// Find which mesh to use for mesh selection
ImgData* MorphoViewer::findSelectMesh()
{
  if(processRunning())
    return 0;
  MorphoGraphX* mgx = mainWindow();
  if(!mgx) {
    Information::out << "Could not find main window" << endl;
    return 0;
  }
  if(mgx->activeMesh() == 0) {
    if(stack1->mesh->showMesh() and
       (stack1->mesh->showMeshPoints() or stack1->mesh->showMeshLines()))
      return stack1;
  }
  if(mgx->activeMesh() == 1) {
    if(stack2->mesh->showMesh() and
       (stack2->mesh->showMeshPoints() or stack1->mesh->showMeshLines()))
    return stack2;
  }
  return (0);
}

// Find selected triangle, respect the clipping planes
void MorphoViewer::findSelectTriangle(ImgData* stk, uint x, uint y, std::vector<uint>& vList, 
                                                               int& label, bool useParentLabel)
{
  clipEnable();
  stk->findSelectTriangle(x, y, vList, label, useParentLabel);
  clipDisable();
}

// Get interface variables
void MorphoViewer::getGUIFlags(QMouseEvent* e)
{
  mousePos = e->pos();

  shiftPressed = (e->modifiers() & Qt::ShiftModifier);
  altPressed = (e->modifiers() & Qt::AltModifier);
  controlPressed = (e->modifiers() & Qt::ControlModifier);

  leftButton = (e->button() == Qt::LeftButton or (e->buttons() & Qt::LeftButton));
  rightButton = (e->button() == Qt::RightButton or (e->buttons() & Qt::RightButton));
}

// Get interface variables
void MorphoViewer::getGUIFlags(QKeyEvent* e)
{
  shiftPressed = (e->modifiers() & Qt::ShiftModifier);
  altPressed = (e->modifiers() & Qt::AltModifier);
  controlPressed = (e->modifiers() & Qt::ControlModifier);
}

void MorphoViewer::keyReleaseEvent(QKeyEvent* e)
{
  getGUIFlags(e);
  if(guiAction == STACK_VOXEL_EDIT) {
    checkVoxelCursor(altPressed);
    updateAll();
  }

  // Cannot depend on getting Qt::Key_Alt
  if(guiActionOn and not altPressed) {
    QMouseEvent ev(QEvent::MouseButtonRelease, QPoint(), Qt::LeftButton, Qt::LeftButton, Qt::AltModifier);
    callGuiAction(guiAction, &ev);
    guiActionOn = false;
  }
  QGLWidget::keyReleaseEvent(e);
}

void MorphoViewer::enterEvent(QEvent* e)
{
  checkVoxelCursor(false);
  setFocus(Qt::MouseFocusReason);
  updateAll();
  QGLWidget::enterEvent(e);
}

void MorphoViewer::leaveEvent(QEvent* e)
{
  checkVoxelCursor(false);
  updateAll();
  QGLWidget::leaveEvent(e);
}

// Keyboard actions, return true if handled, false to pass through
void MorphoViewer::keyPressEvent(QKeyEvent* e)
{
  int key = e->key();
  Qt::KeyboardModifiers mod = e->modifiers();
  float fly = FlySpeed;
  bool handled = false;
  getGUIFlags(e);

  // Grab x,y,z directions in correct frame
  Vec x(1.0, 0.0, 0.0), y(0.0, 1.0, 0.0), z(0.0, 0.0, 1.0), revpt(0, 0, 0);
  if(manipulatedFrame() and manipulatedFrame() != camera()->frame()) {
    x = manipulatedFrame()->transformOf(camera()->frame()->inverseTransformOf(x));
    y = manipulatedFrame()->transformOf(camera()->frame()->inverseTransformOf(y));
    z = manipulatedFrame()->transformOf(camera()->frame()->inverseTransformOf(z));
    fly *= -1;
    // Change revolve aroung pt or not?
  }
  if(key == Qt::Key_Alt and guiAction == STACK_VOXEL_EDIT) {
    checkVoxelCursor(true);
    handled = true;
  } else if(mod == Qt::NoModifier) {
    switch(key) {
    case Qt::Key_Minus:
      FlySpeed *= .9;
      FlySpeed = trim(FlySpeed, .001f, float(M_PI));
      handled = true;
      break;
    case Qt::Key_Right:     // rotate left/right
      manipulatedFrame()->rotateAroundPoint(qglviewer::Quaternion(z, fly), revpt);
      handled = true;
      break;
    case Qt::Key_Left:
      manipulatedFrame()->rotateAroundPoint(qglviewer::Quaternion(z, -fly), revpt);
      handled = true;
      break;
    case Qt::Key_Up:     // move in/out, for translations always use camera frame
      manipulatedFrame()->translate(camera()->frame()->inverseTransformOf(Vec(0.0, 0.0, -fly)));
      handled = true;
      break;
    case Qt::Key_Down:
      manipulatedFrame()->translate(camera()->frame()->inverseTransformOf(Vec(0.0, 0.0, fly)));
      handled = true;
      break;
    case Qt::Key_Delete:
    case Qt::Key_Backspace:
      handled = true;
      emit deleteSelection();
      break;
    }
  } else if(mod == Qt::ShiftModifier) {
    switch(key) {
    case Qt::Key_Plus:
      FlySpeed *= 1.1;
      FlySpeed = trim(FlySpeed, .001f, float(M_PI));
      handled = true;
      break;
    case Qt::Key_Right:     // the translations
      manipulatedFrame()->translate(camera()->frame()->inverseTransformOf(Vec(-fly, 0.0, 0.0)));
      handled = true;
      break;
    case Qt::Key_Left:
      manipulatedFrame()->translate(camera()->frame()->inverseTransformOf(Vec(fly, 0.0, 0.0)));
      handled = true;
      break;
    case Qt::Key_Up:
      manipulatedFrame()->translate(camera()->frame()->inverseTransformOf(Vec(0.0, -fly, 0.0)));
      handled = true;
      break;
    case Qt::Key_Down:
      manipulatedFrame()->translate(camera()->frame()->inverseTransformOf(Vec(0.0, fly, 0.0)));
      handled = true;
      break;
    }
  } else if(mod == Qt::ControlModifier) {
    switch(key) {
    case Qt::Key_Right:     // the other rotations
      manipulatedFrame()->rotateAroundPoint(qglviewer::Quaternion(y, fly), revpt);
      handled = true;
      break;
    case Qt::Key_Left:
      manipulatedFrame()->rotateAroundPoint(qglviewer::Quaternion(y, -fly), revpt);
      handled = true;
      break;
    case Qt::Key_Up:
      manipulatedFrame()->rotateAroundPoint(qglviewer::Quaternion(x, fly), revpt);
      handled = true;
      break;
    case Qt::Key_Down:
      manipulatedFrame()->rotateAroundPoint(qglviewer::Quaternion(x, -fly), revpt);
      handled = true;
      break;
    }
  }
  if(handled)
    updateAll();
  else if(!altPressed)
    QGLViewer::keyPressEvent(e);
}

void MorphoViewer::checkVoxelCursor(bool altPressed)
{
  if(!voxelEditCursor and altPressed and guiAction == STACK_VOXEL_EDIT) {
    ImgData* stk = findSelectStack();
    if(stk and stk->stack->work()->isVisible()) {
      voxelEditCursor = true;
      qApp->setOverrideCursor(QCursor(Qt::BlankCursor));
    }
  } else if(voxelEditCursor and !(altPressed and guiAction == STACK_VOXEL_EDIT)) {
    voxelEditCursor = false;
    qApp->restoreOverrideCursor();
  }
}

int MorphoViewer::findLabel(const Store* store, Point3f start, Point3f dir) const
{
  const Stack* stk = store->stack();
  size_t loop_count = 0;
  start = Point3f(stk->getFrame().coordinatesOf(Vec(start)));
  dir = Point3f(stk->getFrame().transformOf(Vec(dir)));
  start = stk->worldToImagef(start);
  dir = stk->worldToImageVectorf(dir);
  Point3f end = start + dir;
  dir = normalized(dir);
  if(DEBUG) {
    Information::out << "findLabel() - stack coordinates" << endl;
    Information::out << "Start = " << start << " - End = " << end << " - dir = " << dir << endl;
  }
  int incx = dir.x() > 0 ? 1 : -1;
  int incy = dir.y() > 0 ? 1 : -1;
  int incz = dir.z() > 0 ? 1 : -1;
  Point3f incr_l(0, 0, 0);
  if(dir.x() > 0)
    incr_l.x() = 1;
  if(dir.y() > 0)
    incr_l.y() = 1;
  if(dir.z() > 0)
    incr_l.z() = 1;
  Point3u grid_size = stk->size();
  if(fabs(dir.x()) > EPSILON) {
    if(incx > 0) {
      start -= dir * (start.x() / dir.x());
      end -= dir * ((end.x() - grid_size.x()) / dir.x());
    } else {
      start -= dir * (start.x() - grid_size.x()) / dir.x();
      end -= dir * (end.x() / dir.x());
    }
    if(fabs(dir.y()) > EPSILON) {
      if(incy > 0) {
        if(start.y() < 0)
          start -= dir * (start.y() / dir.y());
        if(end.y() > grid_size.y())
          end -= dir * ((end.y() - grid_size.y()) / dir.y());
      } else {
        if(start.y() > grid_size.y())
          start -= dir * (start.y() - grid_size.y()) / dir.y();
        if(end.y() <= 0)
          end -= dir * (end.y() / dir.y());
      }
    }
    if(fabs(dir.z()) > EPSILON) {
      if(incz > 0) {
        if(start.z() < 0)
          start -= dir * (start.z() / dir.z());
        if(end.z() > grid_size.z())
          end -= dir * ((end.z() - grid_size.z()) / dir.z());
      } else {
        if(start.z() > grid_size.z())
          start -= dir * (start.z() - grid_size.z()) / dir.z();
        if(end.z() <= 0)
          end -= dir * (end.z() / dir.z());
      }
    }
  } else if(fabs(dir.y()) > EPSILON) {
    if(incy > 0) {
      start -= dir * (start.y() / dir.y());
      end -= dir * ((end.y() - grid_size.y()) / dir.y());
    } else {
      start -= dir * (start.y() - grid_size.y()) / dir.y();
      end -= dir * (end.y() / dir.y());
    }
    if(fabs(dir.z()) > EPSILON) {
      if(incz > 0) {
        if(start.z() < 0)
          start -= dir * (start.z() / dir.z());
        if(end.z() > grid_size.z())
          end -= dir * ((end.z() - grid_size.z()) / dir.z());
      } else {
        if(start.z() > grid_size.z())
          start -= dir * (start.z() - grid_size.z()) / dir.z();
        if(end.z() <= 0)
          end -= dir * (end.z() / dir.z());
      }
    }
  } else {
    if(incz > 0) {
      start -= dir * (start.z() / dir.z());
      end -= dir * ((end.z() - grid_size.z()) / dir.z());
    } else {
      start -= dir * (start.z() - grid_size.z()) / dir.z();
      end -= dir * (end.z() / dir.z());
    }
  }
  // Check if we cross the grid at all
  if((start.x() < 0 and end.x() < 0)or (start.x() > grid_size.x() and end.x() > grid_size.x())
    or (start.y() < 0 and end.y() < 0) or (start.y() > grid_size.y() and end.y() > grid_size.y())
    or (start.z() < 0 and end.z() < 0) or (start.z() > grid_size.z() and end.z() > grid_size.z())) {
    // No label possible
    return 0;
  }
  // Now, compute the intersection with the clipping planes
  Clip* clips[3] = { c1, c2, c3 };
  double smm[16];
  stk->getFrame().getMatrix(smm);
  Matrix4f sm(smm);
  Point4f hstart(stk->imageToWorld(start));
  hstart[3] = 1.f;
  Point4f hend(stk->imageToWorld(end));
  hend[3] = 1.f;
  for(int i = 0; i < 3; ++i) {
    Clip* c = clips[i];
    if(c->enabled()) {
      // Get the transform from clipping plane to stack frame
      double mm[16];
      c->frame().inverse().getMatrix(mm);
      Matrix4f cm(mm);
      Matrix4f m = sm * cm;
      Point4f c0 = m * c->normalFormPos();
      Point4f c1 = m * c->normalFormNeg();
      // Now perform the tests
      float testS0 = hstart * c0;
      float testE0 = hend * c0;
      float testS1 = hstart * c1;
      float testE1 = hend * c1;
      if((testS0 < 0 and testE0 < 0)or (testS1 < 0 and testE1 < 0)) {
        return 0;                                                     // This is completly clipped
      } else if((testS0 * testE0 < 0)or (testS1 * testE1 < 0))        // Partial clip
      {
        float lambda0 = testS0 / (testS0 - testE0);
        float lambda1 = testS1 / (testS1 - testE1);
        if(lambda0 > lambda1) {
          float a = lambda0;
          lambda0 = lambda1;
          lambda1 = a;
        }
        Point4f newStart = (lambda0 > 0) ? (1 - lambda0) * hstart + lambda0 * hend : hstart;
        Point4f newEnd = (lambda1 < 1) ? (1 - lambda1) * hstart + lambda1 * hend : hend;
        hstart = newStart;
        hend = newEnd;
      }
    }
  }

  hstart /= hstart[3];
  hend /= hend[3];

  start = stk->worldToImagef(Point3f(hstart));
  end = stk->worldToImagef(Point3f(hend));

  Point3f cp = min(Point3f(grid_size - 1u), max(Point3f(0, 0, 0), start));
  Point3i cell(map(floorf, cp));
  if(DEBUG) {
    Information::out << "********* Start computeLine ***************" << endl;
    Information::out << "dir = " << dir << endl;
    Information::out << "start = " << start << " => " << cell << endl << endl;
  }
  // bool add_corners = false;
  size_t max_loop = grid_size.x() + grid_size.y() + grid_size.z();
  const HVecUS& data = store->data();
  while((end - start) * dir > EPSILON and loop_count++ < max_loop 
                                        and stk->offset(cell) < stk->storeSize()) {
    // Check label
    ushort label = data[stk->offset(cell)];
    if(label > 0) {
      if(DEBUG)
        Information::out << "Found label: " << label << endl;
      return label;
    }
    Point3f l(cell);
    Point3f dl(l - start + incr_l);
    if(DEBUG)
      Information::out << "l = " << l << endl;
    Point3f r = divide(dl, dir);
    if(fabs(dir.x()) < EPSILON)
      r.x() = HUGE_VAL;
    if(fabs(dir.y()) < EPSILON)
      r.y() = HUGE_VAL;
    if(fabs(dir.z()) < EPSILON)
      r.z() = HUGE_VAL;
    Point3f ar = fabs(r);
    if(DEBUG)
      Information::out << "dl = " << dl << endl;
    if(DEBUG)
      Information::out << "r = " << r << endl;
    if(ar.x() <= ar.y() and ar.x() <= ar.z()) {
      start += dir * r.x();
      cell.x() += incx;
      if(ar.y() == ar.x())
        cell.y() += incy;
      if(ar.z() == ar.x())
        cell.z() += incz;
    } else if(ar.y() <= ar.x() and ar.y() <= ar.z()) {
      start += dir * r.y();
      cell.y() += incy;
      if(ar.z() == ar.y())
        cell.z() += incz;
    } else {
      start += dir * r.z();
      cell.z() += incz;
    }
    if(DEBUG)
      Information::out << "New start point = " << start << " => " << cell << endl << endl;
  }
  if(DEBUG) {
    if(loop_count >= max_loop) {
      Information::err << "Error, too many loop count" << endl;
    }
    Information::out << "######### End computeLine ###############" << endl;
  }
  return 0;
}
Point3f MorphoViewer::pointUnderPixel(const QPoint &pos, bool &found)
{
  QPoint texPos = QPoint(pos.x(), drawHeight - pos.y() - 1);
  float depth = depthTexture[texPos.y() * drawWidth + texPos.x()];
  found = depth < 1.0;
  Vec point(pos.x(), pos.y(), depth);
  point = camera()->unprojectedCoordinatesOf(point);
  return Point3f(point);
}

void MorphoViewer::mouseDoubleClickEvent(QMouseEvent* e)
{
  if(e->button() == Qt::LeftButton and e->buttons() == (Qt::LeftButton | Qt::RightButton)
     and e->modifiers() == Qt::NoModifier) {
    bool found;
    Point3f pup = pointUnderPixel(e->pos(), found);
    if(found) {
      camera()->setRevolveAroundPoint(Vec(pup));
      setVisualHintsMask(1);
      updateAll();
    }
  } else
    QGLViewer::mouseDoubleClickEvent(e);
}

void MorphoViewer::LabelColorSlot()
{
  setLabel(0);
  updateAll();
}

//  Takes care of GUI others
void MorphoViewer::ResetViewSlot()
{
  Matrix4d m = Matrix4d::identity();
  camera()->setFromModelViewMatrix(m.data());
  camera()->setSceneCenter(Vec(0, 0, 0));
  camera()->setSceneRadius(SceneRadius);
  camera()->showEntireScene();
  stack1->getMainFrame().setFromMatrix(m.data());
  stack1->getTransFrame().setFromMatrix(m.data());
  stack2->getMainFrame().setFromMatrix(m.data());
  stack2->getTransFrame().setFromMatrix(m.data());
  c1->frame().setFromMatrix(m.data());
  c2->frame().setFromMatrix(m.data());
  c3->frame().setFromMatrix(m.data());
  cutSurf->cut->frame().setFromMatrix(m.data());
  _camera->resetZoom();
  updateAll();
}

void MorphoViewer::ReloadShaders()
{
  textureSurfShader.invalidate();
  volumeSurfShader1.invalidate();
  volumeSurfShader2.invalidate();
  indexSurfShader.invalidate();
  raycastingShader1.invalidate();
  raycastingShader2.invalidate();
  updateAll();
}

void MorphoViewer::UpdateLabels()
{
  stack1->LabelColorsChanged = true;
  stack2->LabelColorsChanged = true;
  updateAll();
}

void MorphoViewer::updateAll()
{
//  updateGL();
  update();
}

void MorphoViewer::loadFile(const QString& pth, bool loadStack2, bool loadWork)
{
  QFileInfo fi(pth);
  QString ext = fi.suffix();
  ImgData* stack = (loadStack2) ? stack2 : stack1;

  if(ext == "fct") {
    stack->setColorMap(pth, loadWork);
    updateAll();
  }
}

bool MorphoViewer::saveImageSnapshot(const QString& fileName, QSize finalSize, double oversampling, bool expand,
                                     bool hasGUI)
{
  return saveImageSnapshot(fileName, finalSize, oversampling, expand, hasGUI, 0);
}

bool MorphoViewer::saveImageSnapshot(const QString& fileName, QSize finalSize, double oversampling, bool,
                                     bool hasGUI, QString* error)
{
  if(!finalSize.isValid())
    finalSize = QSize(this->width(), this->height());

  double sizeFac = std::max((double)finalSize.width() / (double)this->width(), (double)finalSize.height() / (double)this->height());
  finalSize = QSize(this->width() * sizeFac, this->height() * sizeFac);
  QSize textureSize = oversampling * finalSize;
  double totalScaleFac = oversampling * sizeFac;
  if(totalScaleFac != 1.0) {
    ImgData::scaleBar.scaleDrawing(totalScaleFac);
    ImgData::colorBar.scaleDrawing(totalScaleFac);
  }
  int savedWidth = drawWidth;
  int savedHeight = drawHeight;
  drawWidth = textureSize.width();
  drawHeight = textureSize.height();
  QGLFramebufferObjectFormat fboFormat;
  fboFormat.setAttachment(QGLFramebufferObject::Depth);
  QGLFramebufferObject fbo(textureSize, fboFormat);
  if(!fbo.bind()) {
    if(hasGUI)
      QMessageBox::critical(0, "Error creating FBO",
                            QString("Unable to create a FBO of size %1x%2. Choose a lower resolution or decrease oversampling.").arg(drawWidth).arg(drawHeight));
    if(error)
      *error = QString("Unable to create a FBO of size %1x%2").arg(drawWidth).arg(drawHeight);
    resizeGL(savedWidth, savedHeight);
    ImgData::scaleBar.restoreScale();
    ImgData::colorBar.restoreScale();
    return false;
  }
  resizeGL(textureSize.width(), textureSize.height());
  baseFboId = fbo.handle();

  preDraw();
  draw(&fbo);
  postDraw();

  baseFboId = 0;
  resizeGL(savedWidth, savedHeight);

  if(!fbo.release()) {
    if(hasGUI)
      QMessageBox::critical(0, "Error releasing FBO",
                            QString("Unable to release the FBO of size %1x%2. Choose a lower resolution or decrease oversampling.").arg(drawWidth).arg(drawHeight));
    if(error)
      *error = QString("Unable to release the FBO of size %1x%2").arg(drawWidth).arg(drawHeight);
    ImgData::scaleBar.restoreScale();
    ImgData::colorBar.restoreScale();
    return false;
  }

  QImage img = fbo.toImage();

  if(totalScaleFac != 1.0) {
    img = img.scaled(finalSize, Qt::KeepAspectRatio, Qt::SmoothTransformation);
    ImgData::scaleBar.restoreScale();
    ImgData::colorBar.restoreScale();
  }

  img.save(fileName);

  return true;
}

static const QString extToFilter[8][2] = { 
  { QString("jpeg"), QString("JPEG") }, { QString("jpg"), QString("JPEG") },
  { QString("png"), QString("PNG") }, { QString("bmp"), QString("BMP") },
  { QString("ppm"), QString("PPM") }, { QString("eps"), QString("EPS") },
  { QString("ps"), QString("PS") }, { QString("xfig"), QString("XFIG") } };

static const QString filters = QObject::tr("JPEG images") + " (*.jpeg *.jpg)" + ";;" 
  + QObject::tr("PNG images") + " (*.png)" + ";;" + QObject::tr("PPM images") 
  + " (*.ppm)" + ";;" + QObject::tr("BMP images") + " (*.bmp)";

void MorphoViewer::recordMovie(bool on)
{
  if(on) {
    QString filename, filetype;
    QString selectedFilter;
    filename = QFileDialog::getSaveFileName(this, "Record movie", "", filters, &selectedFilter);
    if(!filename.isEmpty()) {
      // Find out the type of the selected file
      QString filetype;
      for(int i = 0; i < 8; ++i) {
        QString ext = "." + extToFilter[i][0];
        if(filename.endsWith(ext, Qt::CaseInsensitive)) {
          filetype = extToFilter[i][1];
          filename = filename.left(filename.size() - filetype.size());
          break;
        }
      }
      if(filetype.isEmpty()) {
        int index = selectedFilter.indexOf(' ');
        filetype = selectedFilter.left(index);
      }
      setSnapshotCounter(0);
      setSnapshotFileName(filename);
      setSnapshotFormat(filetype);
      connect(this, SIGNAL(drawFinished(bool)), this, SLOT(saveScreenshot(bool)));
    } else {
      emit recordingMovie(false);
    }
  } else
    disconnect(this, SIGNAL(drawFinished(bool)), this, SLOT(saveScreenshot(bool)));
}

void MorphoViewer::setSceneBoundingBox(const Point3f& bbox)
{
  if(norm(bbox) > 0) {
    SceneRadius = norm(bbox);
    emit changeSceneRadius(SceneRadius);
  }
  updateAll();
}

void MorphoViewer::ClipEnableSlot(Clip* c, bool _val)
{
  if(_val)
    c->enable();
  else
    c->disable();
  updateAll();
}

void MorphoViewer::ClipGridSlot(Clip* c, bool _val)
{
  if(_val)
    c->showGrid();
  else
    c->hideGrid();
  updateAll();
}

void MorphoViewer::ClipWidthSlot(Clip* c, int _val)
{
  c->setWidth(0.5 * SceneRadius * (float)exp(double(_val) / 2000.0));
  if(c->grid() or c->enabled())
    updateAll();
}

// Repeat is used from the label dialog to select a specific color.
// RSS this routine is used only by the label dialog
void MorphoViewer::selectMeshLabel(int label, int repeat, bool replace)
{
  ImgData* stk = findSelectSurf();
  if(stk->mesh->showMesh()) {
    if(label > 0 or repeat > 0) {
      if(replace)
        stk->clearMeshSelect();
      if(stk->mesh->useParents()) {
        stk->selectParent(label, repeat);
        Information::setStatus(QString("Selected parent %1.").arg(label));
      } else {
        stk->selectLabel(label, repeat);
        Information::setStatus(QString("Selected label %1.").arg(label));
      }
    }
  } else
    Information::setStatus("Selection operations require a visible mesh");
}

/*
 * These are the user interaction routines.
 * It would be nice to create a plug-in mechanism for this, so do not change the uniform 
 * interface if possible
 */

// Do pixel editing operations
bool MorphoViewer::stackVoxelEdit(QMouseEvent *e)
{
  // Only use left button
  if(!leftButton)
    return false;

  // Find which surface selected
  ImgData* stk = findSelectStack();
  if(stk == 0)
    return false;

  // First process release
  if(e->type() == QEvent::MouseButtonRelease) {
    stk->voxelEditStop();
    return true;
  }
 
  // Start the editing process
  if(e->type() == QEvent::MouseButtonPress)
    stk->voxelEditStart(clip1, clip2, clip3);

  // Then edit current postition
  if(e->type() == QEvent::MouseButtonPress or e->type() == QEvent::MouseMove) {
    Point3f p(0, 0, 0), x(1, 0, 0), y(0, 1, 0), z(0, 0, 1);
    // Find select point
    bool surf = false;
    if(stk->stack->work()->isVisible() and cutSurf->cut->isVisible()) {
      if(!stk->findSeedPoint(e->pos().x(), e->pos().y(), *cutSurf, p))
        return false;
      else
        surf = true;
    } else {
      p = (Point3f)camera()->unprojectedCoordinatesOf(
                                       Vec(e->pos().x(), e->pos().y(), 0.0), &stk->getFrame());
      x = (Point3f)camera()->unprojectedCoordinatesOf(
                                       Vec(e->pos().x() + 1, e->pos().y(), 0.0), &stk->getFrame());
      y = (Point3f)camera()->unprojectedCoordinatesOf(
                                       Vec(e->pos().x(), e->pos().y() - 1, 0.0), &stk->getFrame());
      z = (Point3f)camera()->unprojectedCoordinatesOf(
                                       Vec(e->pos().x(), e->pos().y(), 1.0), &stk->getFrame());
      x = (x - p).normalize();
      y = (y - p).normalize();
      z = (z - p).normalize();
    }
    pixelRadius = ImgData::VoxelEditRadius * camera()->pixelGLRatio(Vec(p));
    stk->voxelEdit(pixelRadius, p, x, y, z, surf, selectedLabel);
  } else
    return false;

  return true;
}

bool MorphoViewer::stackPickLabel(QMouseEvent *e)
{
  // Only use left button
  if(!leftButton)
    return false;

  // Do nothing for release
  if(e->type() == QEvent::MouseButtonRelease)
    return true;

  MorphoGraphX* mgx = mainWindow();
  int ac = mgx->activeStack();
  if(ac > -1) {
    Stack* stk = mgx->globalProcess()->stack(ac);
    Store* str = stk->currentStore();
    if(str->isVisible() and str->labels()) {
      Vec vorig, vdir;
      camera()->convertClickToLine(e->pos(), vorig, vdir);
      int l = findLabel(str, Point3f(vorig), Point3f(vdir));
      Information::setStatus(QString("Picked stack label %1.").arg(l));
      setLabel(l);
      return true;
    }
  }
  Information::setStatus("You need the active stack to be labeled.");
  return false;
}

bool MorphoViewer::stackDeleteLabel(QMouseEvent *e)
{
  // Only use left button
  if(!leftButton)
    return false;

  // Do nothing for release
  if(e->type() == QEvent::MouseButtonRelease or e->type() == QEvent::MouseMove)
    return true;

  MorphoGraphX* mgx = mainWindow();
  int ac = mgx->activeStack();
  if(ac > -1) {
    Stack* stk = mgx->globalProcess()->stack(ac);
    Store* str = stk->currentStore();
    if(str->isVisible() and str->labels()) {
      Vec vorig, vdir;
      camera()->convertClickToLine(e->pos(), vorig, vdir);
      int l = findLabel(str, Point3f(vorig), Point3f(vdir));
      if(l) {
        QStringList parms;
        parms << QString::number(l) << "0";
        mgx->LaunchProcess("Stack/Segmentation/Fill Label", parms, Process::PROCESS_RUN, false, false);
      }
      return true;
    }
  }
  Information::setStatus("You need the active stack to be labeled.");
  return false;
}

bool MorphoViewer::stackPickFillLabel(QMouseEvent *e)
{
  // Only use left button
  if(!leftButton)
    return false;

  // Do nothing for release
  if(e->type() == QEvent::MouseButtonRelease or e->type() == QEvent::MouseMove)
    return true;

  bool result = false;
  if(stackPickFill) {
    result = stackPickLabel(e);
    if(result) {
      stackPickFill = false;
      return result;
    }
  } else
    result = stackFillLabel(e);

  stackPickFill = true;
  return result;
}

bool MorphoViewer::stackFillLabel(QMouseEvent *e)
{
  // Only use left button
  if(!leftButton)
    return false;

  // Do nothing for release or move
  if(e->type() == QEvent::MouseButtonRelease or e->type() == QEvent::MouseMove)
    return true;

  MorphoGraphX* mgx = mainWindow();
  int ac = mgx->activeStack();
  if(ac > -1) {
    Stack* stk = mgx->globalProcess()->stack(ac);
    Store* str = stk->currentStore();
    if(str->isVisible() and str->labels()) {
      Vec vorig, vdir;
      camera()->convertClickToLine(e->pos(), vorig, vdir);
      int l = findLabel(str, Point3f(vorig), Point3f(vdir));
      if(l) {
        QStringList parms;
        parms << QString::number(l) << QString::number(selectedLabel);
        mgx->LaunchProcess("Stack/Segmentation/Fill Label", parms, Process::PROCESS_RUN, false, false);
        return true;
      }
    }
  }
  Information::setStatus("You need the active stack to be labeled.");
  return false;
}

// Move points when right button is pressed
bool MorphoViewer::meshMovePoints(QMouseEvent *e)
{
  bool cutSurfActive = cutSurf->cut->drawGrid();
  ImgData* stk = findSelectMesh();
  if(!stk)
    return false;

  // First process release
  if(e->type() == QEvent::MouseButtonRelease) {
    stk->updateTriPos(); // Update point positions
  } else if(e->type() == QEvent::MouseButtonPress) {
    // Save old position
    oldPos = e->pos();
    if(cutSurfActive) {
      // Save reference point
      oldVPos = Point3f(camera()->unprojectedCoordinatesOf(
                           Vec(Point3f(oldPos.x(), oldPos.y(), 0)), cutSurf->getFrame()));
    } else if(stk) {   // Process mesh point select
      // Save reference point
      oldVPos = Point3f(camera()->unprojectedCoordinatesOf(
                           Vec(Point3f(oldPos.x(), oldPos.y(), 0)), &stk->getFrame()));
    } else
      return false;
  } else if(e->type() == QEvent::MouseButtonPress or e->type() == QEvent::MouseMove) {
    if(rightButton) {
      // handle the moving of points
      Point3f s(0,0,0);

      // If we hit the control key while moving it will cancel
      if(controlPressed) {
        s.x() = oldPos.x();
        s.y() = oldPos.y();
      } else {
        s.x() = e->pos().x();
        s.y() = e->pos().y();
      }

      Point3f newVPos;
      if(cutSurfActive) {
        // Move bezier points
        newVPos = Point3f(camera()->unprojectedCoordinatesOf(Vec(s), cutSurf->getFrame()));
        forall(const uint u, cutSurf->selectV)
          cutSurf->cut->bezier().bezierV()[u] += Point3d(newVPos - oldVPos);
        cutSurf->cut->hasChanged();
      } else {
        // Move mesh points
        newVPos = Point3f(camera()->unprojectedCoordinatesOf(Vec(s), &stk->getFrame()));
        Point3d deltaP(newVPos - oldVPos);
        #pragma omp parallel for
        for(uint i = 0; i < stk->selectV.size(); i++) {
          vertex v = stk->selectV[i];
          v->pos += deltaP;
        }
        stk->updatePos();
      }
      oldVPos = newVPos;
    }
  }
  return true;
}

// Start mesh selecting/editing operations
bool MorphoViewer::meshSelectPoints(QMouseEvent *e)
{
  bool cutSurfActive = cutSurf->cut->drawGrid();
  ImgData* stk = findSelectMesh();
  if(!stk)
    return false;

  // First process release
  if(e->type() == QEvent::MouseButtonRelease) {
    if(cutSurfActive) {
      // Handle cutting surface
      if(leftButton) {
        for(uint i = 0; i < cutSurf->cut->bezier().bezPoints().x() * cutSurf->cut->bezier().bezPoints().y(); i++) {
          Point3f pos = Point3f(camera()->projectedCoordinatesOf(
                                 Vec(cutSurf->cut->bezier().bezierV()[i]), cutSurf->getFrame()));
          if(selectRect.contains(QPoint(int(pos.x() + .5), int(pos.y() + .5))))
            cutSurf->selectV.insert(i);
        }
      }
    } else if (stk) {
      // Handle mesh, add points to selected set
 
      // Respect clip planes, first find the position of the origin and the normal vector
      Point3f orig1 = Point3f(c1->frame().inverseCoordinatesOf(Vec(0, 0, 0)));
      Point3f orig2 = Point3f(c2->frame().inverseCoordinatesOf(Vec(0, 0, 0)));
      Point3f orig3 = Point3f(c3->frame().inverseCoordinatesOf(Vec(0, 0, 0)));
      // Then the normal
      Point3f norm1 = Point3f(c1->frame().inverseCoordinatesOf(Vec(c1->normal()))) - orig1;
      Point3f norm2 = Point3f(c2->frame().inverseCoordinatesOf(Vec(c2->normal()))) - orig2;
      Point3f norm3 = Point3f(c3->frame().inverseCoordinatesOf(Vec(c3->normal()))) - orig3;

      if(leftButton) {
        // Handle selection on VV mesh
        VtxVec vList;
        const vvGraph &S = stk->mesh->graph();
        if(!S.empty()) {
          forall(const vertex &v, S) {
            Point3f vpos = multiply(Point3f(v->pos), stk->stack->scale());
            Point3f cvpos = Point3f(stk->getFrame().inverseCoordinatesOf(Vec(vpos)));
    
            // Check if clipped, 
            if((!c1->enabled() or (fabs((cvpos - orig1) * norm1) < c1->width()))
               and (!c2->enabled() or (fabs((cvpos - orig2) * norm2) < c2->width()))
               and (!c3->enabled() or (fabs((cvpos - orig3) * norm3) < c3->width()))) {
              Point3f pos = Point3f(camera()->projectedCoordinatesOf(Vec(vpos), &stk->getFrame()));
              if(selectRect.contains(QPoint(int(pos.x() + .5), int(pos.y() + .5))))
                vList.push_back(v);
            }
          }
          if(controlPressed)
            stk->removeSelect(vList);
          else
            stk->addSelect(vList);
        }
      }
    } else
      return false;

    selectRect.setTopRight(oldPos);
    selectRect.setBottomLeft(oldPos);

    return true;
  }
  if(e->type() == QEvent::MouseButtonPress) {
    // Save position and reset selection rectangle
    oldPos = e->pos();
    selectRect.setBottomRight(e->pos());
    selectRect.setTopLeft(e->pos());
    // Clear previous selection if required
    if(cutSurfActive) {
      // Save reference point
      oldVPos = Point3f(camera()->unprojectedCoordinatesOf(
                           Vec(Point3f(oldPos.x(), oldPos.y(), 0)), cutSurf->getFrame()));

      // If shift not pressed clear selection
      if(leftButton and !shiftPressed and !controlPressed)
        cutSurf->clearSelect();
    } else if(stk) {   // Process mesh point select
      // Save reference point
      oldVPos = Point3f(camera()->unprojectedCoordinatesOf(
                           Vec(Point3f(oldPos.x(), oldPos.y(), 0)), &stk->getFrame()));
      // If shift not pressed clear selection
      if(leftButton and !shiftPressed and !controlPressed) {
        stk->clearMeshSelect();
      }
    } else
      return false;
  }
  if(e->type() == QEvent::MouseButtonPress or e->type() == QEvent::MouseMove) {
    if(leftButton) {
      // Handle sizing of select rectangle
      selectRect.setBottomRight(oldPos);
      selectRect.setTopLeft(e->pos());
      selectRect = selectRect.normalized();
    } else 
      return false;
  } else
    return false;

  return true;
}

bool MorphoViewer::meshSelectLasso(QMouseEvent *e)
{
  bool cutSurfActive = cutSurf->cut->drawGrid();

  ImgData* stk = findSelectMesh();
  if(!stk)
    return false;

  // First process release
  if(e->type() == QEvent::MouseButtonRelease) {
    if(cutSurfActive) {
      // Handle cutting surface
      if(leftButton) {
        for(uint i = 0; i < cutSurf->cut->bezier().bezPoints().x() * cutSurf->cut->bezier().bezPoints().y(); i++) {
          Point3f pos = Point3f(camera()->projectedCoordinatesOf(
                                 Vec(cutSurf->cut->bezier().bezierV()[i]), cutSurf->getFrame()));
          if(selectLasso.LassoContains(pos))
            cutSurf->selectV.insert(i);
        }
      }
    } else if (stk) {
      // Handle mesh, add points to selected set
 
      // Respect clip planes, first find the position of the origin and the normal vector
      Point3f orig1 = Point3f(c1->frame().inverseCoordinatesOf(Vec(0, 0, 0)));
      Point3f orig2 = Point3f(c2->frame().inverseCoordinatesOf(Vec(0, 0, 0)));
      Point3f orig3 = Point3f(c3->frame().inverseCoordinatesOf(Vec(0, 0, 0)));
      // Then the normal
      Point3f norm1 = Point3f(c1->frame().inverseCoordinatesOf(Vec(c1->normal()))) - orig1;
      Point3f norm2 = Point3f(c2->frame().inverseCoordinatesOf(Vec(c2->normal()))) - orig2;
      Point3f norm3 = Point3f(c3->frame().inverseCoordinatesOf(Vec(c3->normal()))) - orig3;

      VtxVec vList;
      if(leftButton) {
        // Handle select
        const vvGraph &S = stk->mesh->graph();
        if(!S.empty()) {
          forall(const vertex &v, S) {
            Point3f vpos = multiply(Point3f(v->pos), stk->stack->scale());
            Point3f cvpos = Point3f(stk->getFrame().inverseCoordinatesOf(Vec(vpos)));
  
            // Check if clipped, 
            if((!c1->enabled() or (fabs((cvpos - orig1) * norm1) < c1->width()))
               and (!c2->enabled() or (fabs((cvpos - orig2) * norm2) < c2->width()))
               and (!c3->enabled() or (fabs((cvpos - orig3) * norm3) < c3->width()))) {
              Point3f pos = Point3f(camera()->projectedCoordinatesOf(Vec(vpos), &stk->getFrame()));
              if(selectLasso.LassoContains(pos))
                vList.push_back(v);
            }
          }
          if(controlPressed)
            stk->removeSelect(vList);
          else
            stk->addSelect(vList);
        }
      } 
    } else
      return false;

    selectLasso.Clear();
    selectLasso.AddPoint(Point3f(oldPos.x(),oldPos.y(),0));

    return true;
  }
  if(e->type() == QEvent::MouseButtonPress) {
    // Save position and reset selection rectangle
    oldPos = e->pos();
    selectLasso.Clear();
    selectLasso.AddPoint(Point3f(e->pos().x(),e->pos().y(),0));  
    // Clear previous selection if required
    if(cutSurfActive) {
      // Save reference point
      oldVPos = Point3f(camera()->unprojectedCoordinatesOf(
                           Vec(Point3f(oldPos.x(), oldPos.y(), 0)), cutSurf->getFrame()));

      // If shift not pressed clear selection
      if(leftButton and !shiftPressed and !controlPressed)
        cutSurf->clearSelect();
    } else if(stk) {   // Process mesh point select
      // Save reference point
      oldVPos = Point3f(camera()->unprojectedCoordinatesOf(
                           Vec(Point3f(oldPos.x(), oldPos.y(), 0)), &stk->getFrame()));
      // If shift not pressed clear selection
      if(leftButton and !shiftPressed and !controlPressed) {
        stk->clearMeshSelect();
      }
    } else
      return false;
  }
  if(e->type() == QEvent::MouseButtonPress or e->type() == QEvent::MouseMove) {
    if(leftButton) {
      selectLasso.AddPoint(Point3f(e->pos().x(),e->pos().y(),0));
    } else 
      return false;
  } else
    return false;

  return true;
}

// Pick the label color from the mesh
bool MorphoViewer::meshPickLabel(QMouseEvent *e)
{
  // Only use left button
  if(!leftButton)
    return false;

  // Do nothing for release
  if(e->type() == QEvent::MouseButtonRelease)
    return true;

  // Find selected surface and triangle
  ImgData *stk = findSelectSurf();
  if(stk == 0)
    return false;
  int label = 0;

  std::vector<uint> vlist;
  findSelectTriangle(stk, e->pos().x(), e->pos().y(), vlist, label);
  if(vlist.size() == 0)
    return false;

  // Get labeling flags
  bool showHeat = stk->mesh->surfView() == Mesh::SURF_LABEL and stk->mesh->surfLabelView() == "Label Heat";

  if(!showHeat and label >= 0) {
    setLabel(label);
    Information::setStatus(QString("Picked %1 %2").arg(stk->mesh->useParents() ? "parent " : "label ").arg(label));
  } else if(showHeat and label > 0) {
    if(stk->mesh->labelHeat().count(label) > 0)
      Information::setStatus(QString("Label %1 has heat value: %2.").arg(label).arg(stk->mesh->labelHeat()[label]));
  } else
    return false;

  return true;
}

// Add a new seed for each click
bool MorphoViewer::meshAddSeed(QMouseEvent *e)
{
  // Only use left button
  if(!leftButton)
    return false;

  // Do nothing for release
  if(e->type() == QEvent::MouseButtonRelease)
    return true;

  // Find selected surface and triangle
  ImgData *stk = findSelectSurf();
  if(stk == 0)
    return false;
  int label = 0;
  std::vector<uint> vlist;
  findSelectTriangle(stk, e->pos().x(), e->pos().y(), vlist, label);
  if(vlist.size() == 0)
    return false;

  // Do not seed over an existing label
  if(label > 0)
    return false;

  // If shift pressed use current label
  if(shiftPressed or e->type() ==  QEvent::MouseMove)
    label = stk->mesh->viewLabel();
  else
    label = stk->mesh->nextLabel();

  // Add to parent table if required
  if(stk->mesh->useParents())
    stk->mesh->parents()[label] = label;

  stk->addSeed(label, vlist);
  setLabel(label);

  if(stk->mesh->useParents())
    Information::setStatus(QString("Label and parent %1 added.").arg(label));
  else
    Information::setStatus(QString("Label %1 added.").arg(label));

  return true;
}

// Add the current seed
bool MorphoViewer::meshCurrentSeed(QMouseEvent *e)
{
  // Only use left button
  if(!leftButton)
    return false;

  // Do nothing for release
  if(e->type() == QEvent::MouseButtonRelease)
    return true;

  // Find selected surface and triangle
  ImgData *stk = findSelectSurf();
  if(stk == 0)
    return false;
  int label = 0;
  std::vector<uint> vlist;
  findSelectTriangle(stk, e->pos().x(), e->pos().y(), vlist, label);
  if(vlist.size() == 0)
    return false;

  //if(stk->mesh->useParents()) // Soeren: commented out to avoid labeling cell and parent when only parent is intended
  //  stk->mesh->parents()[selectedLabel] = selectedLabel;
  stk->addSeed(selectedLabel, vlist);
  if(stk->mesh->useParents())
    Information::setStatus(QString("Label and parent %1 seeded.").arg(selectedLabel));
  else
    Information::setStatus(QString("Label %1 seeded.").arg(selectedLabel));

  return true;
}

// Draw the signal on the mesh
bool MorphoViewer::meshDrawSignal(QMouseEvent *e)
{
  // Only use left button
  if(!leftButton)
    return false;

  // Do nothing for release
  if(e->type() == QEvent::MouseButtonRelease)
    return true;

  // Find selected surface and triangle
  ImgData *stk = findSelectSurf();
  if(stk == 0)
    return false;
  int label = 0;
  std::vector<uint> vlist;
  findSelectTriangle(stk, e->pos().x(), e->pos().y(), vlist, label);
  if(vlist.size() == 0)
    return false;

  bool found = false;
  Point3f pup = pointUnderPixel(e->pos(), found);
  if(!found) {
    Information::out << "Point under pixel not found" << endl;
    return false;
  }

  // If control is pressed we'll erase
  controlPressed = (e->modifiers() & Qt::ControlModifier);
  shiftPressed = (e->modifiers() & Qt::ShiftModifier);
  int signal = controlPressed ? 0 : 65535;

  // If shift pressed do all 3, else find closest vertex
  if(shiftPressed)
    stk->drawSignal(signal, vlist);
  else {
    uint v;
    double distance = std::numeric_limits<double>::max();
    for(auto u : vlist) {
      vertex w = vertex(stk->idVA[u]);
      Point3f wpos = multiply(Point3f(w->pos), stk->stack->scale());
      Point3f pos = Point3f(stk->getFrame().inverseCoordinatesOf(Vec(wpos)));
      float d = norm(pos - pup);
      if(distance > d) {
        distance = d;
        v = u;
      }
    }
    stk->drawSignal(signal, {v});
  }

  return true;
}


// Fill selected label with current label
bool MorphoViewer::meshFillLabel(QMouseEvent *e)
{
  static IntSet labels;

  // Only use left button
  if(!leftButton)
    return false;

  // Find selected surface and triangle
  ImgData *stk = findSelectSurf();
  if(stk == 0)
    return false;

  // Update on release
  if(e->type() == QEvent::MouseButtonRelease) {
    if(labels.size() > 0) {
      if(stk->mesh->useParents())
        stk->setParent(labels, selectedLabel);
      else
        stk->fillLabel(labels, selectedLabel);
    }
    return true;
  } 

  // We'll keep labels in a set to make selection fast
  if(e->type() == QEvent::MouseButtonPress)
    labels.clear();

  int label = 0;
  std::vector<uint> vlist;
  findSelectTriangle(stk, e->pos().x(), e->pos().y(), vlist, label);
  if(vlist.size() == 0)
    return false;

  if(label <= 0)
    return false;

  if(stk->mesh->useParents()) {
    // If filling with parents on we need the actual label
    findSelectTriangle(stk, e->pos().x(), e->pos().y(), vlist, label, false);
    if(labels.count(label) == 0) {
      Information::setStatus(QString("Area with label %1 set to parent %2.").arg(label).arg(selectedLabel));
      labels.insert(label);
    }
  } else {
    if(label != selectedLabel and labels.count(label) == 0) {
      Information::setStatus(QString("Area with label %1 filled with label %2.").arg(label).arg(selectedLabel));
      labels.insert(label);
    }
  }

  return true;
}

// Pick and fill labels 
bool MorphoViewer::meshPickFillLabel(QMouseEvent *e)
{
  // Only use left button
  if(!leftButton)
    return false;

  // Find selected surface and triangle
  ImgData *stk = findSelectSurf();
  if(stk == 0)
    return false;

  if(e->type() == QEvent::MouseButtonPress) {
    int label = 0;
    std::vector<uint> vlist;
    findSelectTriangle(stk, e->pos().x(), e->pos().y(), vlist, label);
    if(vlist.size() == 0 or label <= 0) {
      meshPickFill = true;
      return false;
    }
    if(meshPickFill) {
      setLabel(label);
      meshPickFill = false;
    } else {
      meshPickFill = true;
      IntSet labels;
      labels.insert(label);
      if(stk->mesh->useParents()) {
        // If filling with parents on we need the actual label
        findSelectTriangle(stk, e->pos().x(), e->pos().y(), vlist, label, false);
        if(label <= 0)
          return false;
        stk->setParent(labels, selectedLabel);
      } else
        stk->fillLabel(labels, selectedLabel);
    }
  }
  return true;
}


// Pick and fill labels 
bool MorphoViewer::meshPickLabelFillParent(QMouseEvent *e)
{

  // Do nothing for release
  if(e->type() == QEvent::MouseButtonRelease)
    return true;

  MorphoGraphX* mgx = mainWindow();
  if(!mgx) {
    Information::out << "Could not find main window" << endl;
    return 0;
  }

  if(leftButton){
    // pick label
    meshPickLabel(e);

    // change active stack
    mgx->toggleStackTabs();
  } else {
    // set parent
    ImgData *stk = findSelectSurf();
    if(stk == 0)
      return false;
    stk->mesh->setUseParents(true);

    int label = 0;
    std::vector<uint> vlist;
    findSelectTriangle(stk, e->pos().x(), e->pos().y(), vlist, label);
    if(vlist.size() == 0)
      return false;

    stk->addSeed(selectedLabel, vlist);
  }

  return true;
}

// Add the current label to the selection
bool MorphoViewer::meshSelectLabel(QMouseEvent *e)
{
  static IntSet labels;

  // Only use left button
  if(!leftButton)
    return false;

  // Find selected surface and triangle
  ImgData *stk = findSelectSurf();
  if(stk == 0)
    return false;

  // Update on release
  if(e->type() == QEvent::MouseButtonRelease) {
    if(labels.size() > 0) {
      if(controlPressed) {
        if(stk->mesh->useParents())
          stk->unselectParent(labels);
        else
          stk->unselectLabel(labels);
      } else {
        if(stk->mesh->useParents())
          stk->selectParent(labels);
        else
          stk->selectLabel(labels);
      }
      return true;
    } else 
      return false;
  }
  // We'll keep labels in a set to make selection fast
  if(e->type() == QEvent::MouseButtonPress) {
    labels.clear();
    // Clear selection if shift not pressed
    if(!shiftPressed) {
      stk->clearMeshSelect();
    }
  }

  int label = 0;
  std::vector<uint> vlist;
  findSelectTriangle(stk, e->pos().x(), e->pos().y(), vlist, label);
  if(vlist.size() == 0)
    return false;

  if(label <= 0)
    return false;

  // Insert label and print message
  if(labels.count(label) == 0) {
    if(stk->mesh->useParents()) 
      Information::setStatus(QString("Selected parent %1.").arg(label));
    else
      Information::setStatus(QString("Selected label %1.").arg(label));

    labels.insert(label);
    setLabel(label);
  }

  return true;
}

// Add a connected region to the selection
bool MorphoViewer::meshSelectConnected(QMouseEvent *e)
{
  // Only use left button
  if(!leftButton)
    return false;

  // Do nothing for release
  if(e->type() == QEvent::MouseButtonRelease)
    return true;

  // Find selected surface and triangle, print error msg only on press
  ImgData *stk = findSelectSurf();
  if(stk == 0)
    return false;
  int label = 0;
  std::vector<uint> vlist;
  findSelectTriangle(stk, e->pos().x(), e->pos().y(), vlist, label);
  if(vlist.size() == 0)
    return false;

  // Select connected areas
  if(!shiftPressed and !controlPressed)
    stk->clearMeshSelect();
  stk->selectConnected(vlist, controlPressed);

  return true;
}

// Add a triangle to the selection
bool MorphoViewer::meshSelectTriangle(QMouseEvent *e)
{
  // Only use left button
  if(!leftButton or e->type() == QEvent::MouseButtonRelease)
    return false;

  // Grab the stack
  ImgData *stk = findSelectSurf();
  if(stk == 0)
    return false;

  // Clear existing selection if required
  if(e->type() == QEvent::MouseButtonPress) {
    if(!shiftPressed and !controlPressed) {
      stk->clearMeshSelect();
    }
  }
  if(e->type() == QEvent::MouseButtonPress or e->type() == QEvent::MouseMove) {
    int label = 0;
    std::vector<uint> vlist;
    findSelectTriangle(stk, e->pos().x(), e->pos().y(), vlist, label);
    if(vlist.size() == 0)
      return false;
    stk->selectVertices(vlist, controlPressed);
  }

  return true;
}

// Grab the label from the other mesh to use as seed
bool MorphoViewer::meshGrabSeed(QMouseEvent *e)
{
  // Only use left button
  if(!leftButton)
    return false;

  // Save the label from MousePress for MouseMove, on Release it is cleared
  static int olabel = 0;
  if(e->type() == QEvent::MouseButtonRelease) {
    olabel = 0;
    return true;
  }

  // Find selected surface and triangle, print error msg only on press
  ImgData *stk = findSelectSurf();
  if(stk == 0)
    return false;
  int label = 0;
  std::vector<uint> vlist;
  findSelectTriangle(stk, e->pos().x(), e->pos().y(), vlist, label);
  if(vlist.size() == 0)
    return false;

  // Get other stack
  ImgData* ostk = stack1;
  if(stk == stack1)
    ostk = stack2;
  if(ostk->mesh->empty())
    return false;
  // Find center of selected triangle
  // RSS is center used? 
  Point3d center(0.0, 0.0, 0.0);
  forall(uint u, vlist) {
    vertex v = stk->idVA[u];
    center += v->pos;
  }
  center /= float(vlist.size());

  if(e->type() == QEvent::MouseButtonPress) {
    // Find selected triangle in select frame buffer
    std::vector<uint> ovlist;
    findSelectTriangle(ostk, e->pos().x(), e->pos().y(), ovlist, olabel);
    if(ovlist.size() == 0)
      return false;
  }

  // Return if unlabeled in other mesh
  if(olabel <= 0) {
    setLabel(0);
    return false;
  }

  // If grabbing parent add new seed (and map) or just map if already colored
  if(stk->mesh->useParents()) {
    // re-grab the actual label rather than the parent
    findSelectTriangle(stk, e->pos().x(), e->pos().y(), vlist, label, false);
    if(label > 0)
      stk->setParent(label, olabel);
  } else {
    // Just grabbing the seed
    stk->addSeed(olabel, vlist);
  }
  setLabel(olabel);

  if(stk->mesh->useParents())
    Information::setStatus(QString("Parent label %1 grabbed.").arg(olabel));
  else if(label > 0)
    Information::setStatus(QString("Seed %1 grabbed.").arg(olabel));

  return true;
}

bool MorphoViewer::callGuiAction(GuiAction guiAction, QMouseEvent *e)
{
  // Right button moves points
  if(rightButton and guiAction != MESH_PICK_FILL_PARENT)
    return meshMovePoints(e);

  // Left button is used for editing operation
  switch(guiAction) {
    case STACK_PICK_LABEL:
      // Select 3D label;
      return stackPickLabel(e);
    case STACK_PICK_FILL_LABEL:
      // Pick and Fill 3D label;
      return stackPickFillLabel(e);
    case STACK_FILL_LABEL:
      // Fill 3D label;
      return stackFillLabel(e);
    case STACK_VOXEL_EDIT:
      // 3D voxel editing with magic wand
      return stackVoxelEdit(e);
    case STACK_DEL_LABEL:
      // Delete 3D label;
      return stackDeleteLabel(e);
    case MESH_SEL_POINTS:
      // Mesh point selection
      return meshSelectPoints(e);
    case MESH_SEL_LASSO:
      // Mesh select lasso
      return meshSelectLasso(e);
    case MESH_PICK_LABEL:
      // Pick the label color
      return meshPickLabel(e);
    case MESH_ADD_SEED:
      // Add a new seed on the mesh
      return meshAddSeed(e);
    case MESH_CURR_SEED:
      // Add the current seed
      return meshCurrentSeed(e);
    case MESH_DRAW_SIGNAL:
      // Draw the signal on the mesh
      return meshDrawSignal(e);
    case MESH_GRAB_SEED:
      // Grab the label from the other mesh and use for a seed 
      return meshGrabSeed(e);
    case MESH_FILL_LABEL:
      // Fill with a label
      return meshFillLabel(e);
    case MESH_PICK_FILL_LABEL:
      // Pick and Fill the label
      return meshPickFillLabel(e);
    case MESH_PICK_FILL_PARENT:
      // Pick and Fill the parent label
      return meshPickLabelFillParent(e);
    case MESH_SEL_LABEL:
      // Select by label
      return meshSelectLabel(e);
    case MESH_SEL_CONN:
      // Select connected regions
      return meshSelectConnected(e);
    case MESH_SEL_TRIS:
      // Select connected regions
      return meshSelectTriangle(e);
  }
  return false;
}

void MorphoViewer::mousePressEvent(QMouseEvent *e)
{
  lastButtons = e->buttons();
  getGUIFlags(e);

  // Alt-left button for select/edit
  bool handled = false;
  if(altPressed) 
    handled = callGuiAction(guiAction, e);

  if(handled) {
    guiActionOn = true;
    updateAll();
  }

  if(!altPressed)
    QGLViewer::mousePressEvent(e);
}

void MorphoViewer::mouseMoveEvent(QMouseEvent* e)
{
  getGUIFlags(e);
  bool handled = false;

  // Check for missed events
  if(guiActionOn and (not altPressed or not (leftButton or rightButton))) {
    QMouseEvent ev(QEvent::MouseButtonRelease, e->pos(), Qt::LeftButton, 0, 0);
    callGuiAction(guiAction, &ev);
    guiActionOn = false;
  }

  // Check 3D editing cursor
  if(guiAction == STACK_VOXEL_EDIT) {
    checkVoxelCursor(altPressed);
    if(altPressed)
      handled = true;
  }

  // Alt left button for select/edit
  if(guiActionOn and altPressed)
    handled = callGuiAction(guiAction, e);

  if(handled)
    updateAll();
  if(!altPressed) 
    QGLViewer::mouseMoveEvent(e);
} 

void MorphoViewer::mouseReleaseEvent(QMouseEvent *e)
{
  lastButtons = e->buttons();
  getGUIFlags(e);

  bool handled = false;
  if(guiActionOn) {
    handled = callGuiAction(guiAction, e);
    guiActionOn = false;
  }

  if(handled) {
    updateAll();
  }

  // Always pass through
  QGLViewer::mouseReleaseEvent(e);
}


