//
// This file is part of MorphoGraphX - http://www.MorphoGraphX.org
// Copyright (C) 2012-2015 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.
// 
//
// Note this file is not used. It stored the processes before the plugin system was created. There are still a few
// here that never got converted.
//
#include "ImageData.hpp"
#include "CutSurf.hpp"
#include "ClipRegion.hpp"
#include "Progress.hpp"
#include "Information.hpp"
#include "Misc.hpp"
#include "Insert.hpp"
#include "Vector.hpp"

using namespace std;

#ifdef WIN32
#  define FileDialogOptions QFileDialog::DontUseNativeDialog
#else
const int FileDialogOptions = 0;
#endif

int ImgData::ProcessCommandIndex;
QTreeWidgetItem* ImgData::ProcessCommandItem;
QTreeWidget* ImgData::ProcessTab;
std::vector<string> ImgData::ProcessCommands;
std::vector<Point12f> ImgData::ProcessParms;
std::vector<std::vector<QString> > ImgData::ProcessLabels;
std::map<QTreeWidgetItem*, int> ImgData::ProcessItemMap;
std::map<int, QTreeWidgetItem*> ImgData::ProcessIndexMap;

// Launch a mesh or stack processing operation
int ImgData::processCommand(uint idx, ImgData& stk, CutSurf& cutSurf, ClipRegion& clip1, ClipRegion& clip2,
                            ClipRegion& clip3)
{
  if(idx > ProcessParms.size())
    return 0;

  HVecUS& data = currentData();
  TransferFunction& currentTransferFct();
  Point12f& procP = ProcessParms[idx];
  QString cmd = ProcessLabels[idx][0];
  QApplication::setOverrideCursor(QCursor(Qt::WaitCursor));
  // Result flag tells what to reload
  uint retval = RELOAD_NONE;
  bool cmdCheck = true;
  // Note: parameters start at procP[2]. procP[0] is stack switch policy, and procP[1] is #parms
  try {
    // Stack menu
    if(cmd == "Average") {     // DONE
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processAverage(data, uint(procP[2]), uint(procP[3]), uint(procP[4]), uint(procP[5]));
    } else if(cmd == "Annihilate") {     // DONE
      if(checkProcess(idx, CHKPRC_STACK_AND_MESH))
        retval |= processAnnihilate(data, procP[2], uint(procP[3]));
    } else if(cmd == "Apply transfer function") {     // DONE
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processApplyTransferFunction(data, procP[2], procP[3], procP[4], procP[5]);
    } else if(cmd == "Average Main and Work Stacks") {     // DONE
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processAverageStacks();
    } else if(cmd == "Blob Detect") {     // DONE
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processBlobDetect(data, false, uint(procP[2]));
    } else if(cmd == "Blob Watershed") {     // DONE
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processBlobDetect(data, true, uint(procP[2]));
    } else if(cmd == "Blur") {     // DONE
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processBlur(data, uint(procP[2]));
    } else if(cmd == "Brighten/Darken") {     // DONE
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processBrighten(data, procP[2]);
    } else if(cmd == "Clear Work Stack") {     // DONE
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processClearWorkStack(procP[2]);
    } else if(cmd == "Clip Stack") {     // DONE
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processClipStack(data, clip1, clip2, clip3);
    } else if(cmd == "Color Gradient") {     // DONE
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processColorGrad(data, procP[2]);
    } else if(cmd == "Consolidate Regions") {
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processConsolidateRegions(procP[2]);
    } else if(cmd == "Copy Main to Work Stack") {     // DONE
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processCopyToWork();
    } else if(cmd == "Copy Work to Main Stack") {     // DONE
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processCopyToMain();
    } else if(cmd == "Dilate") {     // DONE
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processDilate(data, uint(procP[2]), uint(procP[3]), uint(procP[4]));
    } else if(cmd == "Erode") {     // DONE
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processErode(data, uint(procP[2]), uint(procP[3]), uint(procP[4]));
    } else if(cmd == "Edge Detect") {
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processEdgeDetect(data, procP[2], procP[3], procP[4], uint(procP[5]));
    } else if(cmd == "Filter") {     // DONE
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processFilter(data, uint(procP[2]), uint(procP[3]));
    } else if(cmd == "Fill Holes") {
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processFillHoles(data, uint(procP[2]), uint(procP[3]), uint(procP[4]), uint(procP[5]));
    } else if(cmd == "Fill Stack to Mesh") {
      if(checkProcess(idx, CHKPRC_STACK_AND_MESH))
        retval |= processTrimStack(data, true, uint(procP[2]));
    } else if(cmd == "Reaction-Diffusion Activator-Inhibitor") {
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processActInh(data, procP);
    } else if(cmd == "Reaction-Diffusion Turing") {
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processTuring(data, procP);
    } else if(cmd == "Invert") {     // DONE
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processInvert(data);
    } else if(cmd == "Laplace Transform") {
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processLaplace(data);
    } else if(cmd == "Local Maxima") {
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processLocalMaxima(data, uint(procP[2]), uint(procP[3]), uint(procP[4]), uint(procP[5]));
    } else if(cmd == "Normalize Stack") {
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processNormalize(data, uint(procP[2]), uint(procP[3]), uint(procP[4]));
    } else if(cmd == "Relabel pixels") {
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processRelabel(data, uint(procP[2]), uint(procP[3]));
    } else if(cmd == "Remove Outlying Labels") {
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processRemoveOutliers(data, uint(procP[2]), uint(procP[3]), uint(procP[4]));
    } else if(cmd == "Resize Canvas") {
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processResizeCanvas(uint(procP[2]), uint(procP[2]), uint(procP[3]));
    } else if(cmd == "Scale Stack") {
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processScaleStack(uint(procP[2]), uint(procP[2]), uint(procP[3]));
    } else if(cmd == "Shift Stack") {
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processShift(uint(procP[2]), uint(procP[3]), uint(procP[4]));
    } else if(cmd == "Subtract Main and Work Stacks") {
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processSubtractStacks();
    } else if(cmd == "Swap Stacks") {
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processSwapStacks();
    } else if(cmd == "Trim Stack") {
      if(checkProcess(idx, CHKPRC_STACK_AND_MESH))
        retval |= processTrimStack(data, false, 0);
    } else if(cmd == "Watershed") {
      if(checkProcess(idx, CHKPRC_STACK))
        retval |= processWatershed();
      // Mesh Menu
    } else if(cmd == "Balloon Segmentation") {
      if(checkProcess(idx, CHKPRC_STACK_AND_MESH))
        retval |= processBalloonSegmentation(data, procP);
    } else if(cmd == "Convert to Cells") {
      if(checkProcess(idx, CHKPRC_LABEL))
        retval |= processConvertCells(procP[2]);
    } else if(cmd == "Curvature Gaussian") {
      if(checkProcess(idx, CHKPRC_STACK_AND_SURF))
        retval |= processCurvature(uint(procP[2]), procP[3], procP[4], "gaussian");
    } else if(cmd == "Curvature SumSquare") {
      if(checkProcess(idx, CHKPRC_STACK_AND_SURF))
        retval |= processCurvature(uint(procP[2]), procP[3], procP[4], "sumsquare");
    } else if(cmd == "Cutting Surf Mesh") {
      if(checkProcess(idx, CHKPRC_DRAWMESH))
        retval |= processMakeSurf(cutSurf);
    } else if(cmd == "Delete Low Heat Labels") {
      if(checkProcess(idx, CHKPRC_LABEL_OR_HEAT))
        retval |= processDeleteLowHeatLabels(procP[2]);
    } else if(cmd == "Fix Corners") {
      if(checkProcess(idx, CHKPRC_LABEL))
        retval |= processFixCorners();
    } else if(cmd == "Heat Map") {
      if(checkProcess(idx, CHKPRC_LABEL_OR_HEAT))
        retval |= processHeatMap(stk, procP[2]);
    } else if(cmd == "Marching Cubes Surface") {
      if(checkProcess(idx, CHKPRC_STACK_AND_DRAWMESH))
        retval |= processMarchCubes(data, false, procP[2], uint(procP[3]));
    } else if(cmd == "Marching Cubes 3D") {
      if(checkProcess(idx, CHKPRC_STACK_AND_DRAWMESH))
        retval |= processMarchCubes(data, true, procP[2], uint(procP[3]));
    } else if(cmd == "Mesh from Local Maxima") {
      if(checkProcess(idx, CHKPRC_STACK_AND_DRAWMESH))
        retval |= processMeshLocalMaxima(data, procP[2]);
    } else if(cmd == "Project Signal") {
      if(checkProcess(idx, CHKPRC_STACK_AND_SURF))
        retval |= processProject(data, false, procP[2], procP[3], 0, 0);         // Done
    } else if(cmd == "Project Signal Absolute") {
      if(checkProcess(idx, CHKPRC_STACK_AND_SURF))
        retval |= processProject(data, true, procP[2], procP[3], procP[4], procP[5]);         // Done
    } else if(cmd == "Segment Clear") {
      if(checkProcess(idx, CHKPRC_LABEL))
        retval |= processSegmentClear();         // Done
    } else if(cmd == "Segment Mesh") {
      if(checkProcess(idx, CHKPRC_LABEL))
        retval |= processSegmentMesh(uint(procP[2]));         // Done
    } else if(cmd == "Smooth Mesh") {
      if(checkProcess(idx, CHKPRC_MESH))
        retval |= processSmoothMesh();
    } else if(cmd == "Smooth Mesh Signal") {
      if(checkProcess(idx, CHKPRC_SURF))
        retval |= processSmoothSignal(procP[2]);
    } else if(cmd == "Shrink Mesh") {
      if(checkProcess(idx, CHKPRC_MESH))
        retval |= processShrinkMesh(procP[2]);
    } else if(cmd == "Subdivide") {
      if(checkProcess(idx, CHKPRC_MESH))
        retval |= processSubdivide();         // Done
    } else if(cmd == "Subdivide Adaptive") {
      if(checkProcess(idx, CHKPRC_MESH))
        retval |= processSubdivideAdapt(procP[2], procP[3]);
    } else
      cmdCheck = false;
  }
  catch(CImgInstanceException& e) {
    SETSTATUS("Error:CImg");
  }
  catch(string mesg) {
    SETSTATUS("Error:" << mesg.data());
  }
  catch(...) {
    SETSTATUS("process:Error:Unknown exception");
  }
  Progress::stop();
  QApplication::restoreOverrideCursor();

  // Reload changed data
  if(retval & RELOAD_MAIN)
    reloadTex(mainDataTexId);
  if(retval & RELOAD_WORK)
    reloadTex(workDataTexId);

  if(retval & RELOAD_VBO) {
    chkGraph();
    fillVBOs();
  } else {   // VBOs includes labels and positions
    if(retval & RELOAD_TRIS)
      updateTris();
    if(retval & RELOAD_LINES)
      updateLines();
    if(retval & RELOAD_POS)
      updatePos();
  }

  // ProcessdParms[0] = 0-do nothing, 1-go to main stack, 2-go to work stack, 3-turn on heatmap
  if(cmdCheck)
    return (int(procP[0]));
  else
    return (0);
}

// Stack and mesh processing operation checks
bool ImgData::checkProcess(uint idx, CheckProcessEnum e)
{
  if(e == CHKPRC_STACK) {
    if(!valid or !(MainShow or WorkShow)) {
      SETSTATUS("Operation requires Main or Work stack to be loaded and visible");
      return (false);
    }
  } else if(e == CHKPRC_MESH) {
    if(S.empty() or !(MeshShow or SurfShow)) {
      SETSTATUS("Operation requires Mesh or Surface to be loaded and visible");
      return (false);
    }
  } else if(e == CHKPRC_SURF) {
    if(S.empty() or !SurfShow) {
      SETSTATUS("Operation requires Surface to be loaded and visible");
      return (false);
    }
  } else if(e == CHKPRC_DRAWMESH) {
    if(!(MeshShow or SurfShow)) {
      SETSTATUS("Operation requires Mesh or Surface to be visible");
      return (false);
    }
  } else if(e == CHKPRC_LABEL) {
    if(S.empty() or !SurfShow or !SurfLabels) {
      SETSTATUS("Operation requires Surface to be loaded and labels visible");
      return (false);
    }
  } else if(e == CHKPRC_NORMAL) {
    if(S.empty() or !SurfShow or SurfLabels or SurfHeat) {
      SETSTATUS("Operation requires Surface to be loaded and visible");
      return (false);
    }
  } else if(e == CHKPRC_LABEL_OR_HEAT) {
    if(S.empty() or !SurfShow or !(SurfLabels or SurfHeat)) {
      SETSTATUS("Operation requires Surface to be loaded and labels or heat visible");
      return (false);
    }
  } else if(e == CHKPRC_STACK_AND_MESH) {
    if(!valid or !(MainShow or WorkShow) or S.empty() or !(MeshShow or SurfShow)) {
      SETSTATUS("Operation requires Stack to be loaded and visible and Mesh or Surface to be loaded and visible");
      return (false);
    }
  } else if(e == CHKPRC_STACK_AND_SURF) {
    if(!valid or S.empty() or !SurfShow) {
      SETSTATUS("Operation requires Stack to be loaded and Surface to be loaded and visible");
      return (false);
    }
  } else if(e == CHKPRC_STACK_AND_DRAWMESH) {
    if(!valid or !(MainShow or WorkShow) or !(MeshShow or SurfShow)) {
      SETSTATUS("Operation requires Stack to be loaded and visible and Mesh or Surface to be visible");
      return (false);
    }
  }

  uint parms = ProcessParms[idx][1];
  if(parms == 0)
    SETSTATUS(Section << " " << ProcessLabels[idx][0].toStdString());
  else {
    stringstream ss;
    ss << Section << " " << ProcessLabels[idx][0].toStdString() << "(" << ProcessLabels[idx][2].toStdString() << "="
       << ProcessParms[idx][2];
    for(int i = 1; i < parms; i++)
      ss << ", " << ProcessLabels[idx][i + 2].toStdString() << "=" << ProcessParms[idx][i + 2];
    SETSTATUS(ss.str() << ")");
  }

  return (true);
}

#define off(x, y) ((y)*SizeX + (x))
// Fill in holes in an edge detected surface
uint ImgData::processFillHoles(HVecUS& data, uint xradius, uint yradius, uint threshold, uint fillval)
{
  fillval = clip(fillval, uint(0), uint(0xFFFF));

  // First create height map, allocate update map
  std::vector<int> hmap(SizeX * SizeY);
  std::vector<int> umap(SizeX * SizeY);

  for(uint x = 0; x < SizeX; x++) {
    for(uint y = 0; y < SizeY; y++) {
      uint idx = off(x, y);
      hmap[idx] = umap[idx] = 0;
      for(uint z = 0; z < SizeZ; z++)
        if(data[offset(x, y, z)] >= threshold)
          hmap[idx] = z;
    }
  }

  // Scan in X direction
  if(xradius > 0) {
    for(uint y = 0; y < SizeY; y++) {
      // first look for range of interest
      int firstx = 0;
      while(firstx < int(SizeX) and hmap[off(firstx, y)] == 0)
        firstx++;
      if(firstx >= int(SizeX - 1))
        continue;
      int lastx = SizeX - 1;
      while(lastx > firstx and hmap[off(lastx, y)] == 0)
        lastx--;
      if(--lastx <= ++firstx)
        continue;

      for(int x = firstx; x < lastx; x++) {
        uint idx = off(x, y);
        // max prev and next z's in range
        int pz = hmap[idx], nz = hmap[idx];
        int px = x, nx = x;
        int startx = x - xradius;
        int endx = x + xradius;
        if(startx < firstx - 1)
          startx = firstx - 1;
        if(endx > lastx + 1)
          endx = lastx + 1;
        for(int ix = startx; ix < x; ix++)
          if(hmap[off(ix, y)] > pz) {
            pz = hmap[off(ix, y)];
            px = ix;
          }
        for(int ix = endx; ix > x; ix--)
          if(hmap[off(ix, y)] > nz) {
            nz = hmap[off(ix, y)];
            nx = ix;
          }
        // Calculate line intersection with pixel
        if(px < nx) {
          int newz = pz + int(float(x - px) * float(nz - pz) / float(nx - px) + .5);
          if(newz > hmap[idx] and newz > umap[idx]) {
            umap[idx] = newz;
          }
        }
      }
    }
  }

  // Scan in Y direction
  if(yradius > 0) {
    for(uint x = 0; x < SizeX; x++) {
      // first look for range of interest
      int firsty = 0;
      while(firsty < int(SizeY) and hmap[off(x, firsty)] == 0)
        firsty++;
      if(firsty >= int(SizeY - 1))
        continue;
      int lasty = SizeY - 1;
      while(lasty > firsty and hmap[off(x, lasty)] == 0)
        lasty--;
      if(--lasty <= ++firsty)
        continue;

      for(int y = firsty; y < lasty; y++) {
        uint idx = off(x, y);
        // max prev and next z's in range
        int pz = hmap[idx], nz = hmap[idx];
        int py = y, ny = y;
        int starty = y - yradius;
        int endy = y + yradius;
        if(starty < firsty - 1)
          starty = firsty - 1;
        if(endy > lasty + 1)
          endy = lasty + 1;
        for(int iy = starty; iy < y; iy++)
          if(hmap[off(x, iy)] > pz) {
            pz = hmap[off(x, iy)];
            py = iy;
          }
        for(int iy = endy; iy > y; iy--)
          if(hmap[off(x, iy)] > nz) {
            nz = hmap[off(x, iy)];
            ny = iy;
          }
        // Calculate line intersection with pixel
        if(py < ny) {
          int newz = pz + int(float(y - py) * float(nz - pz) / float(ny - py));
          if(newz > hmap[idx] and newz > umap[idx]) {
            umap[idx] = newz;
          }
        }
      }
    }
  }

  // Update stack
  for(int x = 0; x < int(SizeX); x++)
    for(int y = 0; y < int(SizeY); y++)
      for(int z = umap[off(x, y)]; z >= 0; z--)
        WorkData[offset(x, y, z)] = fillval;

  return (RELOAD_WORK);
}

// Adaptively detect edge detect of stack
uint ImgData::processEdgeDetect(HVecUS& data, float edgeThresh, float edgeFactor, float edgeMult, uint fillval)
{
  ushort lowthresh = clip(ushort(float(edgeThresh) * edgeFactor), ushort(0), ushort(0xFFFF));
  ushort fill = clip(fillval, uint(0), uint(0xFFFF));
  Point3i imgSize(SizeX, SizeY, SizeZ);
  edgeDetectCU(imgSize, lowthresh, edgeThresh, edgeMult, fill, data, WorkData);
  return (RELOAD_WORK);
}

// 3D surface evaluation function for marching cubes
int ImgData::eval(Point3f p)
{
  Point3i pi = worldToImage(p);
  int x = pi.x();
  int y = pi.y();
  int z = pi.z();

  if(!boundsOK(x, y, z))
    return (0);

  // Decide which stack to use
  HVecUS& data = *marchData;
  return (data[offset(x, y, z)]);
}

// Create mesh "spot" from local maxima
uint ImgData::processMeshLocalMaxima(HVecUS& data, float radius)
{
  if(radius <= 0)
    return (RELOAD_NONE);

  uint count = 0;
  for(uint i = 0; i < data.size(); i++)
    if(data[i] > 0)
      count++;
  if(count == 0)
    return (RELOAD_NONE);

  try {
    // Clear mesh
    S.clear();

    Progress::start("Finding mesh from local maxima:", count);
    radius = worldToAbstract(radius);
    uint step = 0;
    for(uint z = 0; z < SizeZ; z++)
      for(uint y = 0; y < SizeY; y++)
        for(uint x = 0; x < SizeX; x++) {
          uint label = data[offset(x, y, z)];
          if(label == 0)
            continue;
          if(!Progress::advance(step++))
            throw("Operation cancelled by user");

          Point3f pos = imageToWorld(Point3i(x, y, z));

          // top, bottom, left, right, near(front), far(back)
          vertex t, b, l, r, n, f;
          vertex ntl, nbl, nbr, ntr, ftl, fbl, fbr, ftr;
          S.insert(t);
          S.insert(b);
          S.insert(l);
          S.insert(r);
          S.insert(n);
          S.insert(f);
          S.insert(ntl);
          S.insert(nbl);
          S.insert(nbr);
          S.insert(ntr);
          S.insert(ftl);
          S.insert(fbl);
          S.insert(fbr);
          S.insert(ftr);

          t->label = b->label = l->label = r->label = n->label = f->label = label;
          ntl->label = nbl->label = nbr->label = ntr->label = label;
          ftl->label = fbl->label = fbr->label = ftr->label = label;

          t->pos = pos + Point3f(0.0f, radius, 0.0f);
          b->pos = pos + Point3f(0.0f, -radius, 0.0f);
          l->pos = pos + Point3f(-radius, 0.0f, 0.0f);
          r->pos = pos + Point3f(radius, 0.0f, 0.0f);
          n->pos = pos + Point3f(0.0f, 0.0f, radius);
          f->pos = pos + Point3f(0.0f, 0.0f, -radius);

          float radf = radius * pow(1.0 / 3.0, .5);
          ntl->pos = pos + Point3f(-radf, radf, radf);
          nbl->pos = pos + Point3f(-radf, -radf, radf);
          nbr->pos = pos + Point3f(radf, -radf, radf);
          ntr->pos = pos + Point3f(radf, radf, radf);
          ftl->pos = pos + Point3f(-radf, radf, -radf);
          fbl->pos = pos + Point3f(-radf, -radf, -radf);
          fbr->pos = pos + Point3f(radf, -radf, -radf);
          ftr->pos = pos + Point3f(radf, radf, -radf);

          makeNhbd(S, t, ntl, ntr, ftr, ftl);
          makeNhbd(S, b, nbr, nbl, fbl, fbr);
          makeNhbd(S, l, nbl, ntl, ftl, fbl);
          makeNhbd(S, r, ntr, nbr, fbr, ftr);
          makeNhbd(S, n, ntl, nbl, nbr, ntr);
          makeNhbd(S, f, ftr, fbr, fbl, ftl);

          makeNhbd(S, ntl, ntr, t, ftl, l, nbl, n);
          makeNhbd(S, nbl, nbr, n, ntl, l, fbl, b);
          makeNhbd(S, nbr, ntr, n, nbl, b, fbr, r);
          makeNhbd(S, ntr, nbr, r, ftr, t, ntl, n);
          makeNhbd(S, ftl, ftr, f, fbl, l, ntl, t);
          makeNhbd(S, fbl, ftl, f, fbr, b, nbl, l);
          makeNhbd(S, fbr, fbl, f, ftr, r, nbr, b);
          makeNhbd(S, ftr, ftl, t, ntr, r, fbr, f);
        }
  }
  catch(string mesg) {
    SETSTATUS(mesg.data());
  }
  catch(...) {
    SETSTATUS("processMeshLocalMaxima::Error: unknown exception");
  }

  cells = false;
  imgtex = false;
  return (RELOAD_VBO);
}

// Marching cubes
uint ImgData::processMarchCubes(HVecUS& data, bool cubeLabel, float cubeSize, int pixThresh)
{
  if(cubeLabel)
    SETSTATUS(Section << " Marching cubes by label");
  else
    SETSTATUS(Section << " Marching cubes by threshold");

  typedef std::pair<int, int> IntIntPair;
  typedef std::map<int, int> IntIntMap;
  typedef Vector<2, Point3i> Point3i2;
  typedef std::pair<int, Point3i2> IntPoint3i2Pair;
  typedef std::map<int, Point3i2> IntPoint3i2Map;

  try {
    // Clear mesh
    S.clear();
    fillVBOs();

    marchData = &data;

    // Bounding region
    IntIntMap labcnt;
    IntPoint3i2Map labels;
    Progress::start("Finding object bounds-" + Section, SizeX);
    for(uint x = 0; x < SizeX; x++) {
      if(!Progress::advance(x))
        throw("Marching Cubes Cancelled");

      for(uint y = 0; y < SizeY; y++) {
        for(uint z = 0; z < SizeZ; z++) {
          // Count pixels
          int label;
          if(cubeLabel)
            label = data[offset(x, y, z)];
          else {
            // For normal march only count "inside" pixels
            if(data[offset(x, y, z)] < pixThresh)
              continue;
            label = 0;
          }

          labcnt[label]++;

          // Find bounding box
          if(labcnt[label] == 1)
            labels[label][0] = labels[label][1] = Point3i(x, y, z);
          else {
            if(x < labels[label][0].x())
              labels[label][0].x() = x;
            if(x > labels[label][1].x())
              labels[label][1].x() = x;
            if(y < labels[label][0].y())
              labels[label][0].y() = y;
            if(y > labels[label][1].y())
              labels[label][1].y() = y;
            if(z < labels[label][0].z())
              labels[label][0].z() = z;
            if(z > labels[label][1].z())
              labels[label][1].z() = z;
          }
        }
      }

      // Delete labels less than threshold pixels and label numbers less than 2
      // (0 is unlabeled, inr images use 1)
      if(cubeLabel) {
        forall(const IntIntPair& p, labcnt)
          if(p.second < pixThresh || p.first <= 1)
            labels.erase(p.first);
      }
    }

    // Call marching cubes for each label
    if(cubeLabel)
      Progress::start("Marching cubes-" + Section, labels.size());
    int step = 0;

    forall(const IntPoint3i2Pair& p, labels) {
      // Triangle and vertex vectors
      std::vector<vertex> vertices;
      std::vector<Point3i> triangles;

      if(cubeLabel and !Progress::advance(step++))
        throw(string("Marching cubes cancelled"));
      marchLabel = p.first;
      Point3f bBox[2] = { Point3f(imageToWorld(p.second[0])), Point3f(imageToWorld(Point3i(p.second[1]))) };

      // Polygonizer object
      if(pol)
        delete pol;
      pol = new Polygonizer(this, worldToAbstract(cubeSize), false);

      // Call polygonizer, catch any errors
      pol->march(true, bBox);

      // Vertices to std_vetor
      VCount = pol->no_vertices();
      for(int i = 0; i < VCount; i++) {
        vertex v;
        v->pos = pol->get_vertex(i);

        v->label = marchLabel;
        vertices.push_back(v);
      }

      TCount = pol->no_triangles();
      for(int i = 0; i < TCount; i++)
        triangles.push_back(Point3i(pol->get_triangle(i)));

      // Insert vertices and triangles into mesh
      Mesh::meshFromTriangles(vertices, triangles);

      if(cubeLabel)
        SETSTATUS("Label:" << marchLabel << ", triangles added:" << TCount << ", vertices:" << S.size());
      else
        SETSTATUS("Triangles added:" << TCount << ", vertices:" << S.size());
    }
  }
  catch(string mesg) {
    SETSTATUS(mesg.data());
  }
  catch(...) {
    SETSTATUS("processMarchCubes::Error: unknown exception");
  }

  cells = false;
  imgtex = false;
  return (RELOAD_VBO);
}

// Make a planar mesh surface based on the clipping region
uint ImgData::processMakeSurf(CutSurf& cutSurf)
{
  try {
    std::vector<Point3f> points;
    int uSize = 0, vSize = 0;
    cutSurf.getSurfPoints(this, points, uSize, vSize);
    if(points.empty() or uSize < 2 or vSize < 2)
      return (RELOAD_NONE);

    // Clear the current mesh
    S.clear();

    // Save the vertex ids for connecting
    vertex vtx[uSize][vSize];

    // Create the vertices and fill in the position
    for(int u = 0; u < uSize; u++) {
      for(int v = 0; v < vSize; v++) {
        vertex a;
        S.insert(a);
        vtx[u][v] = a;
        a->pos = points[u * vSize + v];
        a->txpos.x() = clipTex(float(u) / float(uSize - 1));
        a->txpos.y() = clipTex(float(v) / float(vSize - 1));
      }
    }

    // Connect neighborhoods
    for(int u = 0; u < uSize; u++) {
      for(int v = 0; v < vSize; v++) {
        std::vector<vertex> nhbs;
        nhbs.clear();
        nhbs.push_back(vtx[u][v]);
        if(u == 0 and v == vSize - 1) {         // top left corner;
          nhbs.push_back(vtx[u][v - 1]);
          nhbs.push_back(vtx[u + 1][v - 1]);
          nhbs.push_back(vtx[u + 1][v]);
        } else if(u == 0 and v == 0) {         // bottom left corner
          nhbs.push_back(vtx[u + 1][v]);
          nhbs.push_back(vtx[u][v + 1]);
        } else if(u == uSize - 1 and v == 0) {         // bottom right corner
          nhbs.push_back(vtx[u][v + 1]);
          nhbs.push_back(vtx[u - 1][v + 1]);
          nhbs.push_back(vtx[u - 1][v]);
        } else if(u == uSize - 1 and v == vSize - 1) {         // top right corner
          nhbs.push_back(vtx[u - 1][v]);
          nhbs.push_back(vtx[u][v - 1]);
        } else if(v == vSize - 1) {         // top edge
          nhbs.push_back(vtx[u - 1][v]);
          nhbs.push_back(vtx[u][v - 1]);
          nhbs.push_back(vtx[u + 1][v - 1]);
          nhbs.push_back(vtx[u + 1][v]);
        } else if(u == 0) {         // left edge
          nhbs.push_back(vtx[u][v - 1]);
          nhbs.push_back(vtx[u + 1][v - 1]);
          nhbs.push_back(vtx[u + 1][v]);
          nhbs.push_back(vtx[u][v + 1]);
        } else if(v == 0) {         // bottom edge
          nhbs.push_back(vtx[u + 1][v]);
          nhbs.push_back(vtx[u][v + 1]);
          nhbs.push_back(vtx[u - 1][v + 1]);
          nhbs.push_back(vtx[u - 1][v]);
        } else if(u == uSize - 1) {         // right edge
          nhbs.push_back(vtx[u][v + 1]);
          nhbs.push_back(vtx[u - 1][v + 1]);
          nhbs.push_back(vtx[u - 1][v]);
          nhbs.push_back(vtx[u][v - 1]);
        } else {         // Interior vertex
          nhbs.push_back(vtx[u + 1][v - 1]);
          nhbs.push_back(vtx[u + 1][v]);
          nhbs.push_back(vtx[u][v + 1]);
          nhbs.push_back(vtx[u - 1][v + 1]);
          nhbs.push_back(vtx[u - 1][v]);
          nhbs.push_back(vtx[u][v - 1]);
        }
        makeNhbd(nhbs);
      }
    }

    cells = false;
  }
  catch(...) {
    SETSTATUS("processMakeSurf::Error unknow exception");
  }
  SETSTATUS(Section << " Mesh made from cutting surface");
  return (RELOAD_VBO);
}

// Do a step of smoothing
uint ImgData::processSmoothMesh()
{
  // We'll use the normal as a temp position to save space in the vertex struct
  forall(const vertex& v, S) {
    if(MeshSelect and MeshPoints and selectV.size() > 0 and !v->selected)
      continue;
    //    forall(uint u, selectV) {
    //      vertex v = pidVA[u];
    //      v->label = 0;
    //    }
    float wt = 0;
    v->nrml = Point3f(0, 0, 0);
    forall(const vertex& n, S.neighbors(v)) {
      if(!v->margin) {
        vertex m = S.nextTo(v, n);
        float a = triangleArea(v->pos, n->pos, m->pos);
        wt += a;
        v->nrml += Point3f(a * (n->pos + m->pos) / 2.0);
      } else if(n->margin) {
        float a = (n->pos - v->pos).norm();
        v->nrml += Point3f(a * n->pos);
        wt += a;
      }
    }
    if(wt > 0)
      v->nrml /= wt;
    else
      v->nrml = v->pos;
  }
  forall(const vertex& v, S) {
    if(MeshSelect and MeshPoints and selectV.size() > 0 and !v->selected)
      continue;
    v->pos = v->pos * .5 + v->nrml * .5;
  }

  imgtex = false;
  return (RELOAD_POS);
}

// Shrink mesh
uint ImgData::processShrinkMesh(float shrinkAmt)
{
  forall(const vertex& v, S)
    v->pos -= worldToAbstract(shrinkAmt) * v->nrml;
  return (RELOAD_POS);
}

// Subdivide triangle mesh (could be cells)
uint ImgData::processSubdivide()
{
  try {
    vvgraph T;
    T = S;
    Insert<vvgraph> insert;

    // Add vertex in middle of all existing edges
    forall(const vertex& n, T) {
      forall(const vertex& m, T.neighbors(n)) {
        if(m < n)
          continue;
        vertex v = insert(n, m, S);
        v->margin = false;
        v->color = (m->color + n->color) / 2.0;
        v->pos = (m->pos + n->pos) / 2.0;
        v->txpos = clipTex((m->txpos + n->txpos) / 2.0);
        if(m->label == n->label)
          v->label = m->label;
      }
    }
    T = S;
    // Now add new edges
    forall(const vertex& v, T) {
      if(T.valence(v) != 2 or v->margin)
        continue;
      forall(const vertex& n, T.neighbors(v)) {
        vertex m = T.nextTo(v, n);
        vertex nv = T.nextTo(n, v);
        vertex mv = T.prevTo(m, v);
        vertex nm = T.nextTo(nv, n);
        vertex mn = T.prevTo(mv, m);
        if(nm != mn)
          continue;
        S.spliceBefore(nv, nm, mv);
        S.spliceAfter(mv, mn, nv);
      }
    }
  }
  catch(...) {
    SETSTATUS("processSubdivide::Error unknow exception");
  }

  cells = false;
  SETSTATUS(Section << " Subdivide triangles, total vertices:" << S.size());
  return (RELOAD_VBO);
}

// Subdivide all boundary triangles over a threshold area
uint ImgData::processSubdivideAdapt(float cellMaxArea, float borderSize)
{
  try {
    // Find triangles on borders
    markBorder(borderSize);

    // Copy to temp graph
    vvgraph T = S;
    forall(vertex v, T)
      forall(vertex n, T.neighbors(v)) {
        vertex m = T.nextTo(v, n);
        if(!uniqueTri(S, v, n, m))
          continue;

        // Subdivide triangles on border and clear label
        if(v->minb == 0 or m->minb == 0 or n->minb == 0)
          continue;
        v->label = n->label = m->label = 0;
        float area = triangleArea(abstractToWorld(v->pos), abstractToWorld(n->pos), abstractToWorld(m->pos));
        if(area < cellMaxArea)
          continue;

        subdivideBisect(v, n, m);
      }
    SETSTATUS(S.size() << " total vertices");
  }
  catch(...) {
    SETSTATUS("processSubdivideAdapt::Error unknown exception");
  }
  cells = false;
  return (RELOAD_VBO);
}

// Fix artifacts in mesh due to ambiguities in graph watershed
uint ImgData::processFixCorners()
{
  // Do not try to fix cornes on cellular mesh (there should not be any)
  if(cells)
    return (RELOAD_NONE);

  try {
    Progress::start("Fixing border triangles for " + Section, 0);
    int tricount = 0, prevtricount = 1000000;
    while(tricount < prevtricount) {
      prevtricount = tricount;
      std::vector<Point3ul> tris;
      forall(vertex v, S) {
        if(v->label == -1 and v->margin) {
          forall(vertex n, S.neighbors(v)) {
            vertex m = S.nextTo(v, n);
            if(!uniqueTri(S, v, n, m) or !n->margin or n->label != -1 or !m->margin or m->label != -1)
              continue;
            tris.push_back(Point3ul(v.id(), n.id(), m.id()));
          }
        }
      }
      // Delete vertices with valence of <= 2, save neighbors so we can remark the margins
      int vdelcnt = 0;
      std::set<vertex> N;
      forall(const Point3ul& tri, tris) {
        for(int i = 0; i < 2; i++) {
          vertex v(tri[i]);
          if(S.valence(v) <= 2) {
            // Add the neighbor, who might now be on the margin
            forall(const vertex& n, S.neighbors(v))
              N.insert(n);
            vdelcnt++;
            S.erase(v);
            // Remove in neighbor set if there
            if(N.contains(v))
              N.erase(v);
          }
        }
      }
      // Check if neighbors of deleted vertices now on the margin
      markMargin(N, true);

      SETSTATUS("Subdividing region corners for " << Section);
      tris.clear();
      // Subdivide triangles that are in the border
      forall(vertex v, S) {
        if(v->label == -1) {
          forall(vertex n, S.neighbors(v)) {
            vertex m = S.nextTo(v, n);
            if(!uniqueTri(S, v, n, m) or n->label != -1 or m->label != -1)
              continue;
            tris.push_back(Point3ul(v.id(), n.id(), m.id()));
          }
        } else if(v->label == 0) {
          bool allborder = true;
          forall(vertex n, S.neighbors(v))
            if(n->label != -1)
              allborder = false;

          if(allborder) {
            forall(vertex n, S.neighbors(v)) {
              vertex m = S.nextTo(v, n);
              tris.push_back(Point3ul(v.id(), n.id(), m.id()));
            }
          }
        }
      }

      // Subdivide border triangles
      forall(const Point3ul& tri, tris) {
        vertex v(tri.x()), n(tri.y()), m(tri.z());
        v->label = m->label = n->label = 0;
        if(m == S.nextTo(v, n)) {
          tricount++;
          subdivideBisect(v, n, m);
        }
      }
      SETSTATUS(vdelcnt << " margin vertices deleted and " << tricount << " border triangles subdivided for "
                        << Section);
    }
  }
  catch(...) {
    SETSTATUS("processFixCorners::Error unknown exception");
  }
  return (RELOAD_VBO);
}

// Queue types for watershed algorithm
typedef std::pair<float, vertex> cvpair;
typedef std::multimap<float, vertex> cvmap;

// Segment mesh with watershed
uint ImgData::processSegmentMesh(uint maxsteps)
{
  // Start progress bar
  Progress::start("Segmenting mesh-" + Section, S.size());

  cvmap Q;
  do {
    Q.clear();

    // Mark all vertices as not in queue
    forall(const vertex& v, S)
      v->inqueue = false;

    // Put all neighbors of labelled vertices into queue
    forall(const vertex& v, S)
      if(v->label > 0)
        forall(const vertex& n, S.neighbors(v))
          if(n->label == 0)
            if(!n->inqueue) {
              n->inqueue = true;
              Q.insert(cvpair(n->color, n));
            }

    // Process queue
    int steps = maxsteps;
    while(Q.size() > 0 and steps-- > 0) {
      cvmap::iterator cvmax = Q.begin();
      vertex vmax = cvmax->second;

      int label = 0;
      bool foundlabel = false, difflabels = false;
      forall(const vertex& n, S.neighbors(vmax)) {
        if(n->label > 0) {
          if(!foundlabel) {
            label = n->label;
            foundlabel = true;
          } else if(label != n->label)
            difflabels = true;
        } else if(n->label == 0 and !n->inqueue) {
          n->inqueue = true;
          Q.insert(cvpair(n->color, n));
        }
      }
      if(foundlabel) {
        if(difflabels)
          vmax->label = -1;
        else
          vmax->label = label;

        forall(const vertex& n, S.neighbors(vmax))
          if(n->label == 0 and !n->inqueue) {
            n->inqueue = true;
            Q.insert(cvpair(n->color, n));
          }
      }

      vmax->inqueue = false;
      Q.erase(cvmax);
    }

    // Find out how far we have progressed
    int count = 0;
    forall(const vertex& v, S)
      if(v->label == 0)
        count++;
    if(!Progress::advance(S.size() - count))
      break;
    updateTris();
    emit viewerUpdate();
    SETSTATUS("Segmenting " << Section << ", vertices left in queue:" << Q.size());
  } while(Q.size() > 0);

  return (RELOAD_TRIS | RELOAD_LINES);
}

// Clear mesh segmentation
uint ImgData::processSegmentClear()
{
  forall(const vertex& v, S)
    v->label = 0;
  return (RELOAD_TRIS | RELOAD_LINES);
}

// Smooth color
uint ImgData::processSmoothSignal(int passes)
{
  for(int i = 0; i < passes; i++) {
    forall(const vertex& v, S) {
      float tarea = 0;
      v->tcol = 0;
      forall(const vertex& n, S.neighbors(v)) {
        vertex m = S.nextTo(v, n);
        if(!S.edge(m, n))
          continue;
        float a = triangleArea(v->pos, n->pos, m->pos);
        tarea += a;
        v->tcol += a * (v->color + n->color + m->color) / 3.0;
      }
      if(tarea > 0)
        v->tcol /= tarea;
      else
        v->tcol = v->color;
    }
    forall(const vertex& v, S)
      v->color = v->tcol;
  }
  return (RELOAD_TRIS);
}

// Project signal onto mesh
uint ImgData::processProject(HVecUS& data, bool useAbsSignal, float mindist, float maxdist, float absSignalMin,
                             float absSignalMax)
{
  float avgSignal = 0, maxSignal = 0, minSignal = 1e10;
  forall(const vertex& v, S) {
    projSignal(data, v, mindist, maxdist);
    avgSignal += v->signal;
    if(v->signal > maxSignal)
      maxSignal = v->signal;
    if(v->signal < minSignal)
      minSignal = v->signal;
  }
  if(useAbsSignal) {
    minSignal = absSignalMin;
    maxSignal = absSignalMax;
  }
  // Now write color to vertices
  forall(const vertex& v, S)
    v->color = clip((v->signal - minSignal) / (maxSignal - minSignal), 0.0f, 1.0f);

  SETSTATUS(Section << " Signal projection, Avg:" << avgSignal / float(S.size()) << " min:" << minSignal
                    << " max:" << maxSignal);
  return (RELOAD_TRIS);
}

// Get signal and areas by label
void ImgData::getLabelMaps(float borderSize)
{
  int label;

  LabelGeom.clear();
  LabelBordGeom.clear();
  LabelIntGeom.clear();
  LabelTotalSig.clear();
  LabelBordSig.clear();
  LabelIntSig.clear();
  LabelCenter.clear();
  LabelTris.clear();
  LabelPts.clear();
  LabelHeatAmt.clear();

  WallGeom.clear();
  WallBordGeom.clear();
  WallBordSig.clear();
  WallNhbs.clear();
  WallVId.clear();
  WallHeatAmt.clear();

  // Find border vertices for wall maps or area maps that require it
  bool doborder = false;
  if((HeatMapType == HEAT_AREA and (HeatMapRepBord or HeatMapSignal == HEAT_BORD or HeatMapSignal == HEAT_BORD_RATIO
                                    or HeatMapSignal == HEAT_INT or HeatMapSignal == HEAT_INT_RATIO)))
    doborder = true;

  // Find all vertices within a limited distance from the walls
  if(doborder or HeatMapType == HEAT_WALL)
    markBorder(borderSize);

  if(HeatMapType == HEAT_AREA or HeatMapType == HEAT_VOL) {
    // For volumes, delete labels with open meshes
    std::set<int> LabelDel;

    forall(const vertex& v, S) {
      forall(const vertex& n, S.neighbors(v)) {
        vertex m = S.nextTo(v, n);
        if(!uniqueTri(S, v, n, m))
          continue;
        label = getLabel(v, n, m);
        if(label <= 0)
          continue;

        if(HeatMapType == HEAT_AREA) {
          // Total area, signal, and label center
          float area
            = triangleArea(abstractToWorld(v->pos), abstractToWorld(n->pos), abstractToWorld(m->pos));
          LabelGeom[label] += area;
          LabelCenter[label] += area * Point3f(v->pos + n->pos + m->pos) / 3.0;
          LabelTotalSig[label] += area * (v->signal + n->signal + m->signal) / 3.0;

          // Calculate how much of triangle is in border area
          if(doborder) {
            int border = 0;
            if(v->minb != 0)
              border++;
            if(n->minb != 0)
              border++;
            if(m->minb != 0)
              border++;
            // Border area and signal
            float bamt = float(border) / 3.0;
            LabelBordSig[label] += area * bamt * (v->signal + n->signal + m->signal) / 3.0;
            LabelBordGeom[label] += area * bamt;

            // Interior area and signal
            float iamt = float(3 - border) / 3.0;
            LabelIntSig[label] += area * iamt * (v->signal + n->signal + m->signal) / 3.0;
            LabelIntGeom[label] += area * iamt;
          }
        } else if(HeatMapType == HEAT_VOL) {         // Do volume heat map
          // If any vertices of triangle are not the same as the label, then the label
          // does not define a closed volume, add it to deletion list
          if(v->label != label or n->label != label or m->label != label)
            LabelDel.insert(label);

          // Calculate volume enclosed by mesh and save triangles
          LabelGeom[label]
            += signedTetraVolume(abstractToWorld(v->pos), abstractToWorld(n->pos), abstractToWorld(m->pos));

          // Save triangles and points if processing volume signal
          if(HeatMapSignal == HEAT_TOTAL or HeatMapRepSig) {
            LabelTris[label].push_back(Point3ul(v.id(), n.id(), m.id()));
            LabelPts[label].insert(v.id());
            LabelPts[label].insert(n.id());
            LabelPts[label].insert(m.id());
          }
        }
      }
    }

    // Remove labels which define volume calculations on open meshes
    if(HeatMapType == HEAT_VOL) {
      forall(int i, LabelDel) {
        LabelTotalSig.erase(i);
        SETSTATUS("MakeHeatMap::Warning:" << Section << " open mesh for Vol calc, erasing label " << i);
      }

      // Get mesh volume flourescense
      if(HeatMapSignal == HEAT_TOTAL or HeatMapRepSig) {
        HVecUS& data = MainData;         // currentData();

        // Loop over labels
        uint step = 0;
        Progress::start("Finding signal inside mesh " + Section, LabelGeom.size());
        forall(const IntFloatPair& p, LabelGeom) {
          int label = p.first;

          // Process triangle and point lists
          IntIntMap vmap;
          HVec3F pts(LabelPts[label].size());
          int idx = 0;
          forall(const uint u, LabelPts[label]) {
            vertex v(u);
            pts[idx] = v->pos;
            vmap[u] = idx++;
          }
          // Convert triangle index list to new vertex list
          forall(Point3u& tri, LabelTris[label])
            tri = Point3u(vmap[tri.x()], vmap[tri.y()], vmap[tri.z()]);

          forall(Point3u& tri, LabelTris[label])
            if(tri.x() > pts.size() or tri.y() > pts.size() or tri.z() > pts.size())
              cout << "triangle Error" << endl;

          // Get bounding box
          Point3f bBox[2];
          getbBox(pts, bBox);
          Point3i bBoxImg[2] = { worldToImage(bBox[0]), worldToImage(bBox[1]) };
          HVecU mask;

          // call cuda
          Point3i imgSize(SizeX, SizeY, SizeZ);
          if(!Progress::advance(step++))
            throw(string("Operation canceled by user"));

          insideMeshCU(imgSize, zFactor, bBoxImg, step, pts, LabelTris[label], mask);

          // sum results
          forall(const uint idx, mask) {
            if(idx >= SizeXYZ) {
              cout << "getLabelMaps::Error:Invalid idx:" << idx << " size:" << SizeXYZ << endl;
              continue;
            }
            LabelTotalSig[label] += data[idx];
          }

          // If there are pixels with signal, calculate average
          if(mask.size() > 0) {
            LabelTotalSig[label] /= float(mask.size());
            LabelTotalSig[label] *= LabelGeom[label];
          } else
            LabelTotalSig[label] = 0;
        }
      }
    }

    // Fill in heatmap color
    forall(const IntFloatPair& p, LabelGeom) {
      int label = p.first;

      // Calculate label center
      LabelCenter[label] = (LabelGeom[label] == 0 ? Point3f(0, 0, 0) : LabelCenter[label] / LabelGeom[label]);

      float totalSignal, bordSignal, intSignal;

      // Set signal amounts, average over area if required
      totalSignal = LabelTotalSig[label];
      if(doborder) {
        bordSignal = LabelBordSig[label];
        intSignal = LabelIntSig[label];
      }
      if(HeatMapAvg) {
        totalSignal = (LabelGeom[label] == 0 ? 0 : totalSignal / LabelGeom[label]);
        if(doborder) {
          bordSignal = (LabelBordGeom[label] == 0 ? 0 : bordSignal / LabelBordGeom[label]);
          intSignal = (LabelIntGeom[label] == 0 ? 0 : intSignal / LabelIntGeom[label]);
        }
      }

      // Fill in heat depending on user choice
      switch(HeatMapSignal) {
      case HEAT_GEOM:
        LabelHeatAmt[label] = LabelGeom[label];
        break;
      case HEAT_BORD:
        LabelHeatAmt[label] = bordSignal;
        break;
      case HEAT_BORD_RATIO:
        LabelHeatAmt[label] = (totalSignal == 0 ? 0 : bordSignal / totalSignal);
        break;
      case HEAT_INT:
        LabelHeatAmt[label] = intSignal;
        break;
      case HEAT_INT_RATIO:
        LabelHeatAmt[label] = (totalSignal == 0 ? 0 : intSignal / totalSignal);
        break;
      case HEAT_TOTAL:
        LabelHeatAmt[label] = totalSignal;
        break;
      }
    }
  } else if(HeatMapType == HEAT_WALL) {   // Do Wall heatmap
    // Loop through all pairs of cell wall vertices
    forall(const vertex& v, S) {
      if(v->label != -1)
        continue;
      forall(const vertex& n, S.neighbors(v)) {
        if(n->label != -1 or n > v)
          continue;
        IntSet Sv, Sn;
        // Get list of neighbors for each vertex.
        forall(const vertex& s, S.neighbors(v))
          if(s->label != -1)
            Sv.insert(s->label);
        forall(const vertex& s, S.neighbors(n))
          if(s->label != -1)
            Sn.insert(s->label);
        IntVec LabelV;
        forall(int lab, Sv)
          if(Sn.find(lab) != Sn.end())
            LabelV.push_back(lab);
        if(LabelV.size() <= 1)
          continue;
        // Find wall length, add to neighbors array, and add vertices in wall to set
        for(int i = 0; i < LabelV.size() - 1; i++) {
          IntIntPair wall;

          wall = IntIntPair(LabelV[i], LabelV[i + 1]);
          if(wall.first > 0) {
            WallGeom[wall] += (abstractToWorld(v->pos) - abstractToWorld(n->pos)).norm();
            if(HeatMapPolarity)
              LabelGeom[wall.first] += (abstractToWorld(v->pos) - abstractToWorld(n->pos)).norm();
            WallNhbs[wall.first].insert(wall.second);
            WallVId[wall].insert(v.id());
            WallVId[wall].insert(n.id());
          }

          wall = IntIntPair(LabelV[i + 1], LabelV[i]);
          if(wall.first > 0) {
            WallGeom[wall] += (abstractToWorld(v->pos) - abstractToWorld(n->pos)).norm();
            if(HeatMapPolarity)
              LabelGeom[wall.first] += (abstractToWorld(v->pos) - abstractToWorld(n->pos)).norm();
            WallNhbs[wall.first].insert(wall.second);
            WallVId[wall].insert(v.id());
            WallVId[wall].insert(n.id());
          }
        }
      }
    }

    // Find area, and signal on each wall
    if(HeatMapSignal == HEAT_TOTAL or HeatMapRepSig) {
      forall(const vertex& v, S) {
        forall(const vertex& n, S.neighbors(v)) {
          vertex m = S.nextTo(v, n);
          if(!uniqueTri(S, v, n, m))
            continue;

          // Sum signal and area
          IntIntPair wall;
          int label = getLabel(v, n, m);
          if(isBordTriangle(label, v, n, m, wall)) {
            float area
              = triangleArea(abstractToWorld(v->pos), abstractToWorld(n->pos), abstractToWorld(m->pos));
            WallBordGeom[wall] += area;
            WallBordSig[wall] += (v->signal + n->signal + m->signal) / 3.0 * area;
            if(HeatMapPolarity) {
              LabelBordGeom[wall.first] += area;
              LabelBordSig[wall.first] += (v->signal + n->signal + m->signal) / 3.0 * area;
            }
          }
        }
      }
    }

    forall(const IntIntFloatPair& wall, WallGeom) {
      // Fill in heat depending on user choice
      switch(HeatMapSignal) {
      case HEAT_GEOM:
        // Heat based on wall length
        if(HeatMapPolarity and HeatMapCellAvg)
          WallHeatAmt[wall.first]
            = (LabelGeom[wall.first.first] == 0 ? 0 : WallGeom[wall.first] / LabelGeom[wall.first.first]);
        else
          WallHeatAmt[wall.first] = WallGeom[wall.first];
        break;
      case HEAT_TOTAL:
        // Heat based on wall signal
        WallHeatAmt[wall.first] = WallBordSig[wall.first];
        // Divide by area if required
        if(HeatMapAvg)
          WallHeatAmt[wall.first]
            = (WallBordGeom[wall.first] == 0 ? 0 : WallHeatAmt[wall.first] / WallBordGeom[wall.first]);
        // Do relative to cell average
        if(HeatMapPolarity and HeatMapCellAvg) {
          if(HeatMapAvg)
            WallHeatAmt[wall.first]
              = (LabelBordSig[wall.first.first] == 0 or LabelBordGeom[wall.first.first] == 0
                 ? 0
                 : WallHeatAmt[wall.first]
                 / (LabelBordSig[wall.first.first] / LabelBordGeom[wall.first.first]));
          else
            WallHeatAmt[wall.first]
              = (LabelBordSig[wall.first.first] == 0 ? 0 : WallHeatAmt[wall.first]
                 / LabelBordSig[wall.first.first]);
        }
        break;
      }
    }

    // Maps for polarity calc relative to minimum wall
    if(HeatMapPolarity and HeatMapWallMin) {
      IntFloatMap LabelMin;
      // Find min heat per cell
      forall(const IntIntFloatPair& wall, WallHeatAmt)
        if(LabelMin.find(wall.first.first) == LabelMin.end())
          LabelMin[wall.first.first] = wall.second;
        else if(LabelMin[wall.first.first] > wall.second)
          LabelMin[wall.first.first] = wall.second;
      // Divide heat amounts
      forall(const IntIntFloatPair& wall, WallHeatAmt)
        WallHeatAmt[wall.first]
          = (LabelMin[wall.first.first] == 0 ? 0 : WallHeatAmt[wall.first] / LabelMin[wall.first.first]);
    }
  }
}

// Make a heat map
uint ImgData::processHeatMap(ImgData& img, float borderSize)
{
  // Return value, set on throw
  bool ret = true;
  try {
    // Fill dialog with parms (toggle ones that enable/disable other fields)
    heatMapDialog->HeatMapChange->setChecked(!HeatMapChange);
    heatMapDialog->HeatMapChange->setChecked(HeatMapChange);
    heatMapDialog->HeatMapIncr->setChecked(HeatMapIncr);
    heatMapDialog->HeatMapDecr->setChecked(!HeatMapIncr);
    heatMapDialog->HeatMapAvg->setChecked(HeatMapAvg);
    heatMapDialog->HeatMapPolarity->setChecked(HeatMapPolarity);
    heatMapDialog->HeatMapCellAvg->setChecked(HeatMapCellAvg);
    heatMapDialog->HeatMapWallMin->setChecked(HeatMapWallMin);
    heatMapDialog->HeatMapRatio->setChecked(HeatMapRatio);
    heatMapDialog->HeatMapDiff->setChecked(!HeatMapRatio);
    heatMapDialog->HeatMapRange->setChecked(!HeatMapRange);
    heatMapDialog->HeatMapRange->setChecked(HeatMapRange);

    heatMapDialog->HeatMapRepGeom->setChecked(HeatMapRepGeom);
    heatMapDialog->HeatMapRepSig->setChecked(HeatMapRepSig);
    heatMapDialog->HeatMapRepBord->setChecked(HeatMapRepBord);
    heatMapDialog->RangeMin->setValue(HeatMapManRange.x());
    heatMapDialog->RangeMax->setValue(HeatMapManRange.y());

    if(HeatMapType == HEAT_AREA)
      heatMapDialog->HeatMapTypeArea->setChecked(true);
    else if(HeatMapType == HEAT_VOL)
      heatMapDialog->HeatMapTypeVol->setChecked(true);
    else if(HeatMapType == HEAT_WALL)
      heatMapDialog->HeatMapTypeWall->setChecked(true);

    if(HeatMapSignal == HEAT_GEOM)
      heatMapDialog->HeatMapGeom->setChecked(true);
    else if(HeatMapSignal == HEAT_BORD)
      heatMapDialog->HeatMapBord->setChecked(true);
    else if(HeatMapSignal == HEAT_INT)
      heatMapDialog->HeatMapInt->setChecked(true);
    else if(HeatMapSignal == HEAT_BORD_RATIO)
      heatMapDialog->HeatMapBordRatio->setChecked(true);
    else if(HeatMapSignal == HEAT_INT_RATIO)
      heatMapDialog->HeatMapIntRatio->setChecked(true);
    else if(HeatMapSignal == HEAT_TOTAL)
      heatMapDialog->HeatMapTotal->setChecked(true);

    // Set dialog position
    QPoint qp = parent->pos();
    heatMapDlog->move(qp.x() + 100, qp.y() + 100);

    // Run dialog to get parms from user
    QApplication::restoreOverrideCursor();
    if(!heatMapDlog->exec())
      throw(string("Heat map cancelled by user"));
    QApplication::setOverrideCursor(QCursor(Qt::WaitCursor));

    HeatMapChange = heatMapDialog->HeatMapChange->isChecked();
    HeatMapIncr = heatMapDialog->HeatMapIncr->isChecked();
    HeatMapAvg = heatMapDialog->HeatMapAvg->isChecked();
    HeatMapPolarity = heatMapDialog->HeatMapPolarity->isChecked();
    HeatMapCellAvg = heatMapDialog->HeatMapCellAvg->isChecked();
    HeatMapWallMin = heatMapDialog->HeatMapWallMin->isChecked();
    HeatMapRatio = heatMapDialog->HeatMapRatio->isChecked();
    HeatMapRange = heatMapDialog->HeatMapRange->isChecked();
    HeatMapManRange.x() = heatMapDialog->RangeMin->value();
    HeatMapManRange.y() = heatMapDialog->RangeMax->value();
    if(HeatMapChange and img.S.size() <= 0)
      throw(string("Change map requires both meshes loaded"));

    HeatMapRepGeom = heatMapDialog->HeatMapRepGeom->isChecked();
    HeatMapRepSig = heatMapDialog->HeatMapRepSig->isChecked();
    HeatMapRepBord = heatMapDialog->HeatMapRepBord->isChecked();

    // Heat map geometry type
    if(heatMapDialog->HeatMapTypeArea->isChecked())
      HeatMapType = HEAT_AREA;
    else if(heatMapDialog->HeatMapTypeVol->isChecked())
      HeatMapType = HEAT_VOL;
    else if(heatMapDialog->HeatMapTypeWall->isChecked())
      HeatMapType = HEAT_WALL;

    // Cell meshes can only report geometry
    if(cells) {
      SETSTATUS("makHeatMap::Warning: Can only report geometry for cellular mesh");
      heatMapDialog->HeatMapGeom->setChecked(true);
      heatMapDialog->HeatMapRepSig->setChecked(false);
      heatMapDialog->HeatMapRepBord->setChecked(false);
    }

    // Build analysis type string
    sigStr = "";
    if(HeatMapAvg and !heatMapDialog->HeatMapGeom->isChecked())
      sigStr += "Average ";

    // Heat map signal type
    if(heatMapDialog->HeatMapGeom->isChecked()) {
      HeatMapSignal = HEAT_GEOM;
      if(heatMapDialog->HeatMapTypeArea->isChecked())
        sigStr += "Area";
      else if(heatMapDialog->HeatMapTypeVol->isChecked())
        sigStr += "Volume";
      else if(heatMapDialog->HeatMapTypeWall->isChecked())
        sigStr += "Wall";
    } else if(heatMapDialog->HeatMapBord->isChecked()) {
      HeatMapSignal = HEAT_BORD;
      sigStr += "Border Signal";
    } else if(heatMapDialog->HeatMapInt->isChecked()) {
      HeatMapSignal = HEAT_INT;
      sigStr += "Interior Signal";
    } else if(heatMapDialog->HeatMapBordRatio->isChecked()) {
      HeatMapSignal = HEAT_BORD_RATIO;
      sigStr += "Border/Total Signal";
    } else if(heatMapDialog->HeatMapIntRatio->isChecked()) {
      HeatMapSignal = HEAT_INT_RATIO;
      sigStr += "Interior/Total Signal";
    } else if(heatMapDialog->HeatMapTotal->isChecked()) {
      HeatMapSignal = HEAT_TOTAL;
      sigStr += "Total Signal";
    } else
      throw(string("Cannot determine heatmap type"));

    // Get data from the mesh, areas centers, signal
    getLabelMaps(borderSize);

    // If we are doing a change map, get it from the other mesh and take the ratio as the signal
    // Otherwise just copy the info
    LabelHeat.clear();
    WallHeat.clear();
    if(!HeatMapChange) {
      if(HeatMapType == HEAT_AREA or HeatMapType == HEAT_VOL)
        LabelHeat = LabelHeatAmt;
      else if(HeatMapType == HEAT_WALL)
        WallHeat = WallHeatAmt;
    } else {
      if(HeatMapIncr)
        sigStr += " increase";
      else
        sigStr += " decrease";

      // Get values from other mesh
      img.getLabelMaps(borderSize);

      // Change map for area and volume
      if(HeatMapType == HEAT_AREA or HeatMapType == HEAT_VOL) {
        std::vector<int> LabelDel;
        forall(const IntFloatPair& p, LabelHeatAmt) {
          int i = p.first;
          // Calculate change
          if(img.LabelHeatAmt.find(i) != img.LabelHeatAmt.end()) {
            if(HeatMapIncr and HeatMapRatio)
              LabelHeat[i] = (LabelHeatAmt[i] == 0 ? 0 : img.LabelHeatAmt[i] / LabelHeatAmt[i]);
            else if(!HeatMapIncr and HeatMapRatio)
              LabelHeat[i] = (img.LabelHeatAmt[i] == 0 ? 0 : LabelHeatAmt[i] / img.LabelHeatAmt[i]);
            else if(HeatMapIncr and !HeatMapRatio)
              LabelHeat[i] = img.LabelHeatAmt[i] - LabelHeatAmt[i];
            else if(!HeatMapIncr and !HeatMapRatio)
              LabelHeat[i] = LabelHeatAmt[i] - img.LabelHeatAmt[i];
          } else
            LabelDel.push_back(i);
        }
        // Delete labels not in both meshes
        forall(int i, LabelDel) {
          LabelHeat.erase(i);
          SETSTATUS("MakeHeatMap Warning " << Section << " erasing label " << i);
        }
      } else if(HeatMapType == HEAT_WALL) {       // Calculate wall change map
        std::vector<IntIntPair> WallDel;
        forall(const IntIntFloatPair& p, WallHeatAmt) {
          IntIntPair i = p.first;
          // Calculate change
          if(img.WallHeatAmt.find(i) != img.WallHeatAmt.end()) {
            if(HeatMapIncr and HeatMapRatio)
              WallHeat[i] = (WallHeatAmt[i] == 0 ? 0 : img.WallHeatAmt[i] / WallHeatAmt[i]);
            else if(!HeatMapIncr and HeatMapRatio)
              WallHeat[i] = (img.WallHeatAmt[i] == 0 ? 0 : WallHeatAmt[i] / img.WallHeatAmt[i]);
            else if(HeatMapIncr and !HeatMapRatio)
              WallHeat[i] = img.WallHeatAmt[i] - WallHeatAmt[i];
            else if(!HeatMapIncr and !HeatMapRatio)
              WallHeat[i] = WallHeatAmt[i] - img.WallHeatAmt[i];
            else
              SETSTATUS("MakeHeatMap::Error:" << Section << ", unknown ChangeMap type");
          } else
            WallDel.push_back(i);
        }
        // Delete labels not in both meshes
        forall(IntIntPair i, WallDel) {
          WallHeat.erase(i);
          SETSTATUS("MakeHeatMap::Warning:" << Section << " erasing wall " << i.first << "->" << i.second);
        }
      } else
        throw(string("Invalid heat map type"));

      if(HeatMapRatio)
        sigStr += " (ratio)";
      else
        sigStr += " (difference)";
    }

    // Calculate ranges if not using user defined range
    heatMinAmt = HeatMapManRange.x();
    heatMaxAmt = HeatMapManRange.y();
    if(!HeatMapRange) {
      heatMinAmt = 1e10, heatMaxAmt = -1e10;

      // Find max and min
      if(HeatMapType == HEAT_AREA or HeatMapType == HEAT_VOL) {
        forall(const IntFloatPair& p, LabelHeat) {
          float amt = LabelHeat[p.first];
          if(amt < heatMinAmt)
            heatMinAmt = amt;
          if(amt > heatMaxAmt)
            heatMaxAmt = amt;
        }
      } else if(HeatMapType == HEAT_WALL) {
        forall(const IntIntFloatPair& p, WallHeat) {
          float amt = WallHeat[p.first];
          if(amt < heatMinAmt)
            heatMinAmt = amt;
          if(amt > heatMaxAmt)
            heatMaxAmt = amt;
        }
      }
    }

    // Calculate label colors
    float diff = heatMaxAmt - heatMinAmt;
    if(HeatMapType == HEAT_AREA or HeatMapType == HEAT_VOL) {
      forall(const IntFloatPair& p, LabelHeat)
        LabelHeat[p.first] = (diff == 0 ? 0 : clipTex((LabelHeat[p.first] - heatMinAmt) / diff));
    } else if(HeatMapType == HEAT_WALL) {
      forall(const IntIntFloatPair& p, WallHeat)
        WallHeat[p.first] = (diff == 0 ? 0 : clipTex((WallHeat[p.first] - heatMinAmt) / diff));
    }

    // Find units
    QString gunits(""), hunits("");
    switch(HeatMapType) {
    case (HEAT_WALL):
      gunits = UM;
      break;
    case (HEAT_AREA):
      gunits = UM2;
      break;
    case (HEAT_VOL):
      gunits = UM3;
      break;
    }
    if(!HeatMapChange and HeatMapSignal == HEAT_GEOM)
      hunits = gunits;

    heatMapUnit = hunits;

    // Write to file
    if(HeatMapRepGeom or HeatMapRepSig or HeatMapRepBord) {
      QString msg(string(sigStr + " HeatMap results for " + Section).data());
      QString filename = QFileDialog::getSaveFileName(parent, QString("Select file to save " + msg), "",
                                                      QString("CSV files (*.csv)"), 0, FileDialogOptions);

      // Check ending, add suffix if required
      QString suffix(".csv");
      if(filename.right(suffix.size()) != suffix)
        filename.append(suffix);

      QFile file(filename);
      if(!file.open(QIODevice::WriteOnly | QIODevice::Text)) {
        SETSTATUS("makeHeatMap::Error:Cannot open output file:" << filename.toStdString().data());
      } else {
        QTextStream out(&file);

        // Area heat map
        if(HeatMapType == HEAT_AREA) {
          // Write header
          out << "Label, Heat";
          if(HeatMapRepGeom or HeatMapRepSig or HeatMapRepBord) {
            out << ", Area (um^2)";
            if(HeatMapRepBord) {
              out << ", Border Area (um^2)";
              out << ", Interior Area (um^2)";
            }
          }
          if(HeatMapRepSig or HeatMapRepBord) {
            out << ", Total Signal";
            if(HeatMapRepBord) {
              out << ", Border Signal";
              out << ", Interior Signal";
            }
          }
          out << endl;

          // Write data
          forall(const IntFloatPair& p, LabelHeat) {
            int label = p.first;
            out << label << "," << LabelHeat[label];
            if(HeatMapRepGeom or HeatMapRepSig or HeatMapRepBord) {
              out << "," << LabelGeom[label];
              if(HeatMapRepBord) {
                out << "," << LabelBordGeom[label];
                out << "," << LabelIntGeom[label];
              }
            }
            if(HeatMapRepSig or HeatMapRepBord) {
              out << "," << LabelTotalSig[label];
              if(HeatMapRepBord) {
                out << "," << LabelBordSig[label];
                out << "," << LabelIntSig[label];
              }
            }
            out << endl;
          }
        } else if(HeatMapType == HEAT_VOL) {         // Volume heat map
          // Write header
          out << "Label, Heat";
          if(HeatMapRepGeom or HeatMapRepSig) {
            out << ", Volume (um^3)";
            if(HeatMapRepSig)
              out << ", Volume Signal";
          }
          out << endl;

          // Write data
          forall(const IntFloatPair& p, LabelHeat) {
            int label = p.first;
            out << label << "," << LabelHeat[label];
            if(HeatMapRepGeom or HeatMapRepSig) {
              out << "," << LabelGeom[label];
              if(HeatMapRepSig)
                out << "," << LabelTotalSig[label];
            }
            out << endl;
          }
        } else if(HeatMapType == HEAT_WALL) {
          // Write header
          out << "Label, Neighbor, Heat, Wall Length (um)";
          if(HeatMapSignal == HEAT_TOTAL or HeatMapRepSig)
            out << ", Wall Border Area (um^2), Wall Signal";
          out << endl;

          // Write data
          forall(const IntIntFloatPair& p, WallHeat) {
            IntIntPair wall = p.first;
            out << wall.first << ", " << wall.second << "," << WallHeat[wall] << "," << WallGeom[wall];
            if(HeatMapSignal == HEAT_TOTAL or HeatMapRepSig)
              out << "," << WallBordGeom[wall] << "," << WallBordSig[wall];
            out << endl;
          }
        }
        // Footer same for all heat maps
        out << endl << endl << msg.toStdString().data() << endl;

        out << "Main Stack File: " << StackFile.data() << endl;
        if(HeatMapChange)
          out << "Second Stack File: " << img.StackFile.data() << endl;
        out << "Signal range: " << heatMinAmt << "-" << heatMaxAmt << " " << hunits << endl;
        file.close();
      }
    }

    // Print out stats
    int ncells = LabelGeom.size();
    int acells = LabelTotalSig.size();
    double totageom = 0, totgeom = 0;
    if(HeatMapType == HEAT_AREA or HeatMapType == HEAT_VOL) {
      forall(const IntFloatPair& p, LabelGeom) {
        totgeom += LabelGeom[p.first];
        if(LabelTotalSig.find(p.first) != LabelTotalSig.end())
          totageom += LabelGeom[p.first];
      }
      SETSTATUS(Section << " " << sigStr << " analysis, Cells: " << acells << "/" << ncells << ", " << totageom
                        << "/" << totgeom << " " << gunits.toStdString());
    }
    SETSTATUS("Range min: " << heatMinAmt << " max: " << heatMaxAmt << " " << hunits.toStdString());
  }
  catch(string msg) {
    SETSTATUS(Section << "-" << msg << Section);
    ret = false;
  }
  catch(...) {
    SETSTATUS("Error::makeHeatMap:Unknown exception"
              << "-" << Section);
    ret = false;
  }
  return (RELOAD_TRIS);
}

// Process pixel tranformation commands
typedef CImg<float> CImgF;
typedef vector<float> VecF; // Process pixel tranformation commands
typedef CImg<float> CImgF;
typedef vector<float> VecF;

// Average (blur) pixels
uint ImgData::processAverage(HVecUS& data, uint xradius, uint yradius, uint zradius, uint steps)
{
  if(steps <= 0)
    throw(string("Average pixels, steps must be > 0"));
  Point3u imgSize(SizeX, SizeY, SizeZ);
  Point3i radius(xradius, yradius, zradius);
  Progress::start(Section + " Average pixels", steps);
  processPixCU(PROCESS_PIX_AVG, imgSize, radius, data, WorkData);
  if(!Progress::advance(1))
    throw(string("Average Pixels cancelled"));
  for(int i = 1; i < steps; i++) {
    processPixCU(PROCESS_PIX_AVG, imgSize, radius, WorkData, WorkData);
    if(!Progress::advance(i))
      throw(string("Average Pixels cancelled"));
  }
  return (RELOAD_WORK);
}

// Perform pixel erosion (min filter)
uint ImgData::processErode(HVecUS& data, uint xradius, uint yradius, uint zradius)
{
  Point3u imgSize(SizeX, SizeY, SizeZ);
  Point3i radius(xradius, yradius, zradius);
  processPixCU(PROCESS_PIX_ERODE, imgSize, radius, data, WorkData);
  return (RELOAD_WORK);
}

// Perform pixel dilation (max filter)
uint ImgData::processDilate(HVecUS& data, uint xradius, uint yradius, uint zradius)
{
  Point3u imgSize(SizeX, SizeY, SizeZ);
  Point3i radius(xradius, yradius, zradius);
  processPixCU(PROCESS_PIX_DILATE, imgSize, radius, data, WorkData);
  return (RELOAD_WORK);
}

// Filter stack
uint ImgData::processFilter(HVecUS& data, uint lowFilter, uint highFilter)
{
  int highCnt = 0, lowCnt = 0;

  lowFilter = clip(lowFilter, uint(0), uint(0xFFFF));
  highFilter = clip(highFilter, uint(0), uint(0xFFFF));

  Progress::start("High/Low Filter", 100);
  uint steps = 0;
  uint progress = 1;
  for(uint idx = 0; idx < MainData.size(); idx++) {
    if(!steps--) {
      steps = MainData.size() / 100;
      if(!Progress::advance(progress++))
        break;
    }
    if(highFilter > 0 and data[idx] > highFilter) {
      highCnt++;
      WorkData[idx] = highFilter;
    } else if(data[idx] < lowFilter) {
      lowCnt++;
      WorkData[idx] = 0;
    } else
      WorkData[idx] = data[idx];
  }
  SETSTATUS(Section << " Filtered out " << lowCnt << " low pixels, and " << highCnt << " high pixels");
  return (RELOAD_WORK);
}

// Normalize stack
uint ImgData::processNormalize(HVecUS& data, uint xradius, uint yradius, uint zradius)
{
  HVecUS tdata(data.size());
  Point3u imgSize(SizeX, SizeY, SizeZ);
  Point3i radius(xradius, yradius, zradius);
  processPixCU(PROCESS_PIX_DILATE, imgSize, radius, data, tdata);
  for(uint idx = 0; idx < WorkData.size(); idx++)
    if(tdata[idx] > 0)
      WorkData[idx] = uint(data[idx]) * 0xFFFF / uint(tdata[idx]);
    else
      WorkData[idx] = 0;
  return (RELOAD_WORK);
}

// Blur stack
uint ImgData::processBlur(HVecUS& data, uint radius)
{
  if(data != WorkData)
    memcpy(WorkData.data(), data.data(), data.size() * 2);
  CImgUS image(WorkData.data(), SizeX, SizeY, SizeZ, 1, true);
  image.blur(radius);
  return (RELOAD_WORK);
}

// Invert a stack
uint ImgData::processInvert(HVecUS& data)
{
  for(uint i = 0; i < data.size(); i++)
    WorkData[i] = 0xFFFF - data[i];
  return (RELOAD_WORK);
}

// Resize image canvas
uint ImgData::processResizeCanvas(uint xsz, uint ysz, uint zsz)
{
  if(xsz <= 0 or ysz <= 0 or zsz <= 0)
    return (RELOAD_NONE);

  CImgUS mainImage(MainData.data(), SizeX, SizeY, SizeZ, 1, false);
  CImgUS workImage(WorkData.data(), SizeX, SizeY, SizeZ, 1, false);

  setStackSizes(xsz, ysz, zsz, StepXY, StepZ);
  mainImage.resize(SizeX, SizeY, SizeZ, 1, 0, 0, 0, 0, 0, 0);
  workImage.resize(SizeX, SizeY, SizeZ, 1, 0, 0, 0, 0, 0, 0);
  memcpy(MainData.data(), mainImage.data(), MainData.size() * 2);
  memcpy(WorkData.data(), workImage.data(), WorkData.size() * 2);
  return (RELOAD_WORK | RELOAD_MAIN);
}

// Scale a stack
uint ImgData::processScaleStack(uint xsz, uint ysz, uint zsz)
{
  if(xsz <= 0 or ysz <= 0 or zsz <= 0)
    return (RELOAD_NONE);

  CImgUS mainImage(MainData.data(), SizeX, SizeY, SizeZ, 1, false);
  CImgUS workImage(WorkData.data(), SizeX, SizeY, SizeZ, 1, false);
  mainImage.resize(-xsz, -ysz, -zsz, 1, 5);
  workImage.resize(-xsz, -ysz, -zsz, 1, 5);

  setStackSizes(mainImage.width(), mainImage.height(), mainImage.depth(),
                StepXY * float(SizeX) / float(mainImage.width()), StepZ * float(SizeZ) / float(mainImage.depth()));
  memcpy(MainData.data(), mainImage.data(), MainData.size() * 2);
  memcpy(WorkData.data(), workImage.data(), WorkData.size() * 2);
  return (RELOAD_WORK | RELOAD_MAIN);
}

// Shift an image
uint ImgData::processShift(uint xamt, uint yamt, uint zamt)
{
  CImgUS mainImage(MainData.data(), SizeX, SizeY, SizeZ, 1, false);
  mainImage.shift(xamt, yamt, zamt, 0, 0);
  memcpy(MainData.data(), mainImage.data(), MainData.size() * 2);
  CImgUS workImage(WorkData.data(), SizeX, SizeY, SizeZ, 1, false);
  memcpy(WorkData.data(), workImage.data(), WorkData.size() * 2);
  return (RELOAD_WORK | RELOAD_MAIN);
}

// Color Gradient
uint ImgData::processColorGrad(HVecUS& data, float colorGradDiv)
{
  Point3u imgSize(SizeX, SizeY, SizeZ);
  colorGradCU(imgSize, colorGradDiv, data, WorkData);
  return (RELOAD_WORK);
}

// Laplace Transform
uint ImgData::processLaplace(HVecUS& data)
{
  if(data != WorkData)
    memcpy(WorkData.data(), data.data(), data.size() * 2);
  CImgUS image(WorkData.data(), SizeX, SizeY, SizeZ, 1, true);
  image.laplacian();
  return (RELOAD_WORK);
}

// Watershed
uint ImgData::processWatershed()
{
  Progress::start(Section + " Watershed transform", 0);
  if(!Progress::advance(0))
    throw(string("Watershed cancelled"));
  CImgUS image(SizeX, SizeY, SizeZ, 1, 0xFFFF);
  image -= CImgUS(MainData.data(), SizeX, SizeY, SizeZ, 1, true);
  if(!Progress::advance(0))
    throw(string("Watershed cancelled"));
  CImgUS seeds(WorkData.data(), SizeX, SizeY, SizeZ, 1, true);
  seeds.watershed(image);
  Progress::advance(0);
  return (RELOAD_WORK);
}

// Consolidate regions to fix watershed overseeding
uint ImgData::processConsolidateRegions(float tolerance)
{
  // Make a map of light at shared walls
  typedef std::pair<ushort, ushort> UShortUShortPair;
  typedef std::pair<UShortUShortPair, float> UShortUShortFloatPair;
  typedef std::pair<ushort, float> UShortFloatPair;

  Progress::start(Section + " Consolidating regions", 0);
  uint regions, totregions = 0;
  do {
    std::map<UShortUShortPair, float> WallLight;
    std::map<UShortUShortPair, uint> WallPix;
    std::map<ushort, float> CellIntLight;
    std::map<ushort, float> CellIntPix;

    for(int x = 0; x < SizeX; x++) {
      for(int y = 0; y < SizeY; y++)
        for(int z = 0; z < SizeZ; z++) {
          // Find pixel offset
          size_t off = offset(x, y, z);
          uint label = WorkData[off];

          int xs = int(x - 1);
          int xe = int(x + 1);
          if(xs < 0)
            xs = 0;
          if(xe >= SizeX)
            xe = SizeX - 1;

          int ys = int(y - 1);
          int ye = int(y + 1);
          if(ys < 0)
            ys = 0;
          if(ye >= SizeY)
            ye = SizeY - 1;

          int zs = int(z - 1);
          int ze = int(z + 1);
          if(zs < 0)
            zs = 0;
          if(ze >= SizeZ)
            ze = SizeZ - 1;

          bool interior = true;
          for(int nx = xs; nx <= xe; nx++)
            for(int ny = ys; ny <= ye; ny++)
              for(int nz = zs; nz <= ze; nz++) {
                if(nx == x and ny == y and nz == z)
                  continue;
                size_t noff = offset(nx, ny, nz);
                uint nlabel = WorkData[noff];
                if(label == nlabel)
                  continue;
                interior = false;

                UShortUShortPair pr(label, nlabel);
                UShortUShortPair npr(nlabel, label);

                float light = MainData[off] + MainData[noff];
                float zbias = 0.0f;
                if(nz != z) {
                  zbias = 3.0f;
                  if(nx != x)
                    zbias -= 1.0f;
                  if(ny != y)
                    zbias -= 1.0f;
                }
                light += light * zbias / 4.0f;

                WallLight[pr] += light;
                WallLight[npr] += light;
                WallPix[pr]++;
                WallPix[npr]++;
              }
          if(interior) {
            CellIntLight[label] += MainData[off];
            CellIntPix[label]++;
          }
        }
      if(!Progress::advance(0))
        throw(string("Consolidate regions cancelled"));
    }

    if(WallLight.size() == 0)
      throw(string("Error::Consolidate regions: No cells"));

    // Get average light at wall and for whole cell
    std::map<ushort, float> CellWallLight;     // Average light at cell walls
    std::map<ushort, uint> CellWallPix;
    uint totalpix = 0;
    float avglight = 0;
    forall(const UShortUShortFloatPair& pr, WallLight) {
      CellWallLight[pr.first.first] += pr.second;
      CellWallPix[pr.first.first] += WallPix[pr.first];
      avglight += pr.second;
      totalpix += WallPix[pr.first];
    }
    // Check for light in image
    if(totalpix != 0)
      avglight /= totalpix;
    if(avglight == 0)
      throw(string("Error::Consolidate regions: No light in image"));

    // Calculate total light for cell minus the wall
    // Average light at cell wall not including current wall
    std::map<UShortUShortPair, float> CellBordLight;
    std::map<UShortUShortPair, uint> CellBordPix;
    forall(const UShortUShortFloatPair& pr, WallLight) {
      CellBordLight[pr.first] = CellWallLight[pr.first.first] - pr.second;
      CellBordPix[pr.first] = CellWallPix[pr.first.first] - WallPix[pr.first];
    }
    // Divide by Pix counts
    forall(const UShortUShortFloatPair& pr, WallLight) {
      WallLight[pr.first] /= WallPix[pr.first];
      CellBordLight[pr.first] /= CellBordPix[pr.first];
    }
    forall(const UShortFloatPair& pr, CellWallLight)
      CellWallLight[pr.first] /= float(CellWallPix[pr.first]);
    forall(const UShortFloatPair& pr, CellIntLight)
      CellIntLight[pr.first] /= float(CellIntPix[pr.first]);

    // Count neighbors for each cell
    std::map<ushort, uint> Nhbs;
    forall(const UShortUShortFloatPair& pr, WallLight) {
      Nhbs[pr.first.first]++;
      Nhbs[pr.first.second]++;
    }

    // Choose which labels to merge
    std::map<ushort, ushort> MergeMap;
    std::set<ushort> MergeSet;

    // Add cells with low light interfaces, handle when cells encased by a single cell
    forall(const UShortUShortFloatPair& pr, WallLight) {
      if(pr.first.first == 0 or pr.first.second == 0)
        continue;
      float wallLight = pr.second - CellIntLight[pr.first.first];
      float nbwallLight = pr.second - CellIntLight[pr.first.second];
      if(wallLight < 0)
        wallLight = 0;
      if(nbwallLight < 0)
        nbwallLight = 0;

      // Check if wall is unlikely from both directions
      bool chkout, chkin;
      UShortUShortPair rev(pr.first.second, pr.first.first);
      if(Nhbs[pr.first.first] == 1)
        chkout = wallLight / CellWallLight[pr.first.second] < tolerance;
      else
        chkout = wallLight / CellBordLight[pr.first] < tolerance;

      // If there is only one neighbor, CellBordLight will be zero
      if(Nhbs[pr.first.second] == 1)
        chkin = wallLight / CellWallLight[pr.first.first] < tolerance;
      else
        chkin = wallLight / CellBordLight[rev] < tolerance;

      if(chkin and chkout and MergeSet.find(pr.first.first) == MergeSet.end()) {
        MergeMap[pr.first.first] = pr.first.second;
        MergeSet.insert(pr.first.second);
      }
    }

    // Collapse regions
    for(uint i = 0; i < WorkData.size(); i++) {
      if(WorkData[i] > 0 and MergeMap[WorkData[i]] > 0)
        WorkData[i] = MergeMap[WorkData[i]];
      if(i % SizeX * SizeY == 0 and !Progress::advance(0))
        throw(string("Consolidate regions cancelled"));
    }
    regions = MergeSet.size();
    totregions += regions;
    SETSTATUS(Section << " " << totregions << " regions consolidated");
  } while(regions > 0);
  return (RELOAD_WORK);
}

// Average main and work stacks
uint ImgData::processAverageStacks()
{
  for(uint idx = 0; idx < MainData.size(); idx++)
    WorkData[idx] = ushort((uint(MainData[idx]) + uint(WorkData[idx])) / 2);
  return (RELOAD_WORK);
}

// Apply transfer function to stack data
uint ImgData::processApplyTransferFunction(HVecUS& data, float red, float green, float blue, float alpha)
{
  TransferFunction::Colorf minColor, maxColor;
  TransferFunction& fct = currentTransferFct();
  float vmin = HUGE_VAL, vmax = -HUGE_VAL;
  for(int i = 0; i < data.size(); i++) {
    TransferFunction::Colorf c = fct.rgba(data[i] / 65535.0);
    float value = red * c.r() + green * c.g() + blue * c.b() + alpha * c.a();
    if(value < vmin) {
      vmin = value;
      minColor = c;
    }
    if(value > vmax) {
      vmax = value;
      maxColor = c;
    }
    int val = int(floor(value * 65536));
    WorkData[i] = (val > 65535 ? 65535 : val);
  }
  TransferFunction new_fct;
  new_fct.add_rgba_point(0, minColor);
  new_fct.add_rgba_point(1, maxColor);
  workTransferDlg->setTransferFunction(new_fct);
  return (RELOAD_WORK);
}

// Subtract main and work stacks (absolute value)
uint ImgData::processSubtractStacks()
{
  for(uint idx = 0; idx < MainData.size(); idx++)
    WorkData[idx] = ushort(abs(int(MainData[idx]) - int(WorkData[idx])));
  return (RELOAD_WORK);
}

// Increase/decrease brightness
uint ImgData::processBrighten(HVecUS& data, float brightness)
{
  for(uint idx = 0; idx < MainData.size(); idx++) {
    float pix = float(data[idx]) * brightness;
    if(pix < 0)
      pix = 0;
    else if(pix > 0xFFFF)
      pix = 0xFFFF;
    WorkData[idx] = ushort(pix);
  }
  return (RELOAD_WORK);
}

// Copy Main stack to Work stack
uint ImgData::processCopyToWork()
{
  WorkData = MainData;
  return (RELOAD_WORK);
}

// Copy Work stack to Main stack
uint ImgData::processCopyToMain()
{
  MainData = WorkData;
  return (RELOAD_MAIN);
}

// Swap Main and Work stacks
uint ImgData::processSwapStacks()
{
  HVecUS tdata(MainData);
  MainData = WorkData;
  WorkData = tdata;
  return (RELOAD_WORK | RELOAD_MAIN);
}

// Use activator-inhibitor system to find seeds for watershed
uint ImgData::processActInh(HVecUS& data, Point12f& processParms)
{
  HVecUS tdata(data);
  Point3u imgSize(SizeX, SizeY, SizeZ);
  findSeedsRDCU(FIND_SEEDS_ACT_INH, imgSize, processParms, tdata, WorkData);
  return (RELOAD_WORK);
}

// Use Turing-style reaction-diffusion system to find seeds for watershed
uint ImgData::processTuring(HVecUS& data, Point12f& processParms)
{
  HVecUS tdata(data);
  Point3u imgSize(SizeX, SizeY, SizeZ);
  findSeedsRDCU(FIND_SEEDS_TURING, imgSize, processParms, tdata, WorkData);
  return (RELOAD_WORK);
}

typedef std::pair<Point3u, float> PixPair;
bool pixCompare(const PixPair& i, const PixPair& j) {
  return (i.second > j.second);
}
// Detect blobs
uint ImgData::processBlobDetect(HVecUS& data, bool watershed, uint startlabel)
{
  uint label = startlabel;
  HVecUS bdata(data);

  // Sort pixels
  Progress::start(Section + " Blob Detect, sorting pixels", 3);
  if(!Progress::advance(1))
    throw(string("Blob Detect cancelled"));

  std::vector<PixPair> Pix(SizeXYZ);
  uint i = 0;
  for(uint z = 0; z < SizeZ; z++)
    for(uint y = 0; y < SizeY; y++)
      for(uint x = 0; x < SizeX; x++) {
        Pix[i] = PixPair(Point3u(x, y, z), bdata[i]);
        i++;
      }
  if(!Progress::advance(2))
    throw(string("Blob Detect cancelled"));
  std::sort(Pix.begin(), Pix.end(), pixCompare);

  // Prepare result area
  clearData(WorkData);

  Progress::start(Section + "Blob Detect, finding regions", bdata.size());
  for(int i = 0; i < bdata.size(); i++) {
    // Find pixel offset
    const Point3u& p = Pix[i].first;
    size_t off = offset(p.x(), p.y(), p.z());
    if(bdata[off] == 0)
      continue;

    int xs = int(p.x() - 1);
    int xe = int(p.x() + 1);
    if(xs < 0)
      xs = 0;
    if(xe >= SizeX)
      xe = SizeX - 1;

    int ys = int(p.y() - 1);
    int ye = int(p.y() + 1);
    if(ys < 0)
      ys = 0;
    if(ye >= SizeY)
      ye = SizeY - 1;

    int zs = int(p.z() - 1);
    int ze = int(p.z() + 1);
    if(zs < 0)
      zs = 0;
    if(ze >= SizeZ)
      ze = SizeZ - 1;

    bool hashigher = false;
    bool isback = false;
    int labcount = 0, lastlab = 0;
    for(int x = xs; x <= xe; x++)
      for(int y = ys; y <= ye; y++)
        for(int z = zs; z <= ze; z++) {
          if(x == p.x() and y == p.y() and z == p.z())
            continue;
          size_t nboff = offset(x, y, z);
          if(Pix[i].second <= bdata[nboff]) {
            hashigher = true;
            int w = WorkData[nboff];
            if(!watershed and w == 0xFFFF)
              isback = true;
            else if(w > 0 and w != lastlab) {
              labcount++;
              lastlab = w;
            }
          }
        }
    if(!hashigher and labcount == 0)
      WorkData[off] = label++;
    else if(isback)
      WorkData[off] = 0xFFFF;
    else if(labcount == 1)
      WorkData[off] = lastlab;
    else
      WorkData[off] = 0xFFFF;
    if(label == 0xFFFF)
      throw(string("Blob Detect:too many labels"));
    if(i % 100000 == 0 and !Progress::advance(i))
      throw(string("Blob Detect cancelled"));
  }
  // Clear all the background pixels
  for(int i = 0; i < bdata.size(); i++)
    if(WorkData[i] == 0xFFFF)
      WorkData[i] = 0;

  SETSTATUS(Section << " " << label - startlabel << " blobs found");
  return (RELOAD_WORK);
}

// Find Local Maxima
uint ImgData::processLocalMaxima(HVecUS& data, uint xradius, uint yradius, uint zradius, uint startlabel)
{
  Point3i radius(xradius, yradius, zradius);
  uint label = startlabel;

  HVecUS tdata(data.size());
  Point3u imgSize(SizeX, SizeY, SizeZ);
  processPixCU(PROCESS_PIX_DILATE, imgSize, radius, data, tdata);
  uint count = 0;
  for(uint i = 0; i < data.size(); i++) {
    if(tdata[i] == 0)
      continue;
    if(abs(int(tdata[i]) - int(data[i])) == 0) {
      if(startlabel > 0)
        WorkData[i] = label++;
      else
        WorkData[i] = 0x8FFF;
      count++;
    } else
      WorkData[i] = 0;
    if(startlabel > 0 and label == 0xFFFF)
      throw(string("Too many local maxima to label"));
  }
  SETSTATUS(Section << " " << count << " pixels in local maxima");
  return (RELOAD_WORK);
}

typedef std::pair<ushort, uint> UShortUIntPair;
// Remove labels which are too big or too small
uint ImgData::processRemoveOutliers(HVecUS& data, uint minsize, uint maxsize, uint fillval)
{
  // Find label pixel counts
  if(maxsize == 0)
    maxsize = data.size();
  std::map<ushort, uint> PixCount;
  for(uint idx = 0; idx < data.size(); idx++)
    PixCount[data[idx]]++;
  std::set<ushort> LDelete;
  uint count = 0;
  forall(const UShortUIntPair& pr, PixCount)
    if(pr.second < minsize or pr.second > maxsize) {
      count++;
      LDelete.insert(pr.first);
    }
  for(uint idx = 0; idx < data.size(); idx++)
    if(LDelete.find(data[idx]) != LDelete.end())
      WorkData[idx] = fillval;
    else
      WorkData[idx] = data[idx];
  SETSTATUS(Section << " " << count << " labels removed from outlying regions");
  return (RELOAD_WORK);
}

// Relabel pixels after segmentation
uint ImgData::processRelabel(HVecUS& data, uint start, uint step)
{
  std::set<ushort> Labels;
  for(uint idx = 0; idx < data.size(); idx++)
    Labels.insert(data[idx]);
  std::map<ushort, uint> LabMap;
  uint newlab = start;
  forall(const ushort& lab, Labels) {
    if(lab == 0)
      continue;
    LabMap.insert(UShortUIntPair(lab, newlab));
    newlab += step;
    if(newlab > 0xFFFF)
      throw(string("Relabel pixels:Label too large"));
  }
  for(uint idx = 0; idx < data.size(); idx++)
    WorkData[idx] = LabMap[data[idx]];

  SETSTATUS(Section << " " << LabMap.size() << " labels found");
  return (RELOAD_WORK);
}

// Annialate mesh
uint ImgData::processAnnihilate(HVecUS& data, float minDist, float maxDist)
{
  HVecUS tdata(data);

  clearData(WorkData);
  Progress::start("Annihilating " + Section, S.size());
  int cnt = 1;

  // This should be improved to work better on dense meshes, we don't need to look at every point
  forall(const vertex& v, S) {
    if(!Progress::advance(cnt++))
      throw(string("Annihilate cancelled"));
    forall(const vertex& n, S.neighbors(v)) {
      vertex m = S.nextTo(v, n);
      if(!uniqueTri(S, v, n, m))
        continue;

      Point3f nrml = ((v->nrml + n->nrml + m->nrml) / 3.0);
      nrml.normalize();
      nrml /= Scale;
      Point3f vpos = v->pos / Scale;
      Point3f npos = n->pos / Scale;
      Point3f mpos = m->pos / Scale;

      Point3i maxi(0, 0, 0);
      Point3i mini(SizeX, SizeY, SizeZ);
      Point3f pts[2][3];
      Point3i ptsi[2][3];

      pts[0][0] = vpos + nrml * -worldToAbstract(minDist);
      pts[1][0] = vpos + nrml * -worldToAbstract(maxDist);
      pts[0][1] = npos + nrml * -worldToAbstract(minDist);
      pts[1][1] = npos + nrml * -worldToAbstract(maxDist);
      pts[0][2] = mpos + nrml * -worldToAbstract(minDist);
      pts[1][2] = mpos + nrml * -worldToAbstract(maxDist);

      for(int i = 0; i < 2; i++) {
        for(int j = 0; j < 3; j++) {
          ptsi[i][j] = worldToImage(pts[i][j]);
          for(int k = 0; k < 3; k++) {
            if(ptsi[i][j][k] > maxi[k])
              maxi[k] = ptsi[i][j][k];
            if(ptsi[i][j][k] < mini[k])
              mini[k] = ptsi[i][j][k];
          }
        }
      }
      mini.x()--;
      mini.y()--;
      mini.z()--;
      maxi.x()++;
      maxi.y()++;
      maxi.z()++;

      for(int x = mini.x(); x <= maxi.x(); x++) {
        for(int y = mini.y(); y <= maxi.y(); y++) {
          for(int z = mini.z(); z <= maxi.z(); z++) {
            if(!boundsOK(x, y, z))
              continue;
            float s;
            Point3f intp;
            Point3f r0 = imageToWorld(Point3i(x, y, z));
            if(planeLineIntersect(vpos, nrml, r0, r0 + nrml, s, intp) and s >= 0)
              if((intp - r0).norm() <= worldToAbstract(maxDist))
                WorkData[offset(x, y, z)] = tdata[offset(x, y, z)];
          }
        }
      }
    }
  }
  return (RELOAD_WORK);
}

// Clear work stack
uint ImgData::processClearWorkStack(uint fillval)
{
  if(fillval > 0xFFFF)
    fillval = 0xFFFF;
  for(int i = 0; i < WorkData.size(); i++)
    WorkData[i] = fillval;
  return (RELOAD_WORK);
}

// Trim stack to clipping planes
uint ImgData::processClipStack(HVecUS& data, ClipRegion& clip1, ClipRegion& clip2, ClipRegion& clip3)
{
  // Setup matrices and planes for clip test
  getClipTestData(clip1, clip2, clip3);

  // Call cuda to clip
  Point3i imgSize(SizeX, SizeY, SizeZ);
  clipStackCU(imgSize, zFactor, clipDo, Hpn, data, WorkData);
  reloadTex(workDataTexId);

  SETSTATUS(Section << " clipped to intersection of clipping planes");
  return (RELOAD_WORK);
}

// Trim stack to mesh
uint ImgData::processTrimStack(HVecUS& data, bool fill, float fillval)
{
  SETSTATUS("Trimming " << Section << " to mesh");

  // Data for cuda
  HVec3F pts;        // Points for triangles
  HVec3U tris;       // Triangle vertices
  HVecU mask;        // Indices of vertices inside
  Point3f bBox[2];   // Bounding box

  // Get triangle list and bounding box
  getTriangleList(pts, tris);
  getbBox(pts, bBox);
  Point3i bBoxImg[2] = { worldToImage(bBox[0]), worldToImage(bBox[1]) };

  // call cuda to get stack mask
  Point3i imgSize(SizeX, SizeY, SizeZ);
  insideMeshCU(imgSize, zFactor, bBoxImg, 0, pts, tris, mask);

  // write data back to WorkData
  Progress::start("Updating texture for " + Section, 3);
  HVecUS tdata(data);
  Progress::advance(1);
  clearData(WorkData);
  Progress::advance(2);
  uint sz = WorkData.size();
  forall(const uint idx, mask) {
    if(idx >= sz) {
      cout << "trimStack::Error:Invalid idx:" << idx << " size:" << WorkData.size() << endl;
      continue;
    }
    if(fill)
      WorkData[idx] = fillval;
    else
      WorkData[idx] = tdata[idx];
  }
  Progress::advance(3);
  SETSTATUS("Trimmed " << Section << ", " << sz << " pixels kept");
  return (RELOAD_WORK);
}

// Calculate curvature
uint ImgData::processCurvature(uint numnbs, float mincurv, float maxcurv, string curvtype)
{
  if(mincurv >= maxcurv)
    return (RELOAD_NONE);

  bool sumsquare = curvtype == "sumsquare";
  forall(const vertex& v, S) {
    Curvature curv(curvtype);
    std::set<vertex> nbs;
    nbs.insert(v);
    for(int i = 0; i < numnbs; i++) {
      std::set<vertex> tnbs(nbs);
      forall(const vertex& t, tnbs)
        forall(const vertex& n, S.neighbors(t))
          nbs.insert(n);
    }
    std::vector<Point3f> plist, nlist;
    plist.push_back(v->pos);
    nlist.push_back(v->nrml);
    forall(const vertex& n, nbs) {
      if(n == v)
        continue;
      plist.push_back(n->pos);
      nlist.push_back(n->nrml);
    }
    curv.update(plist, nlist);
    v->signal = curv.get_curvature();
    if(sumsquare)
      v->signal = pow(v->signal, .5f);
  }
  forall(const vertex& v, S)
    v->color = clip((v->signal - mincurv) / (maxcurv - mincurv), 0.0f, 1.0f);

  return (RELOAD_TRIS);
}

// Delete all labels with heat less than a threshold
uint ImgData::processDeleteLowHeatLabels(float thresh)
{
  if(LabelHeat.empty())
    return (RELOAD_NONE);
  std::vector<vertex> T;
  forall(const vertex& v, S)
    if(LabelHeat[v->label] < thresh)
      T.push_back(v);
  if(T.empty())
    return (RELOAD_NONE);
  forall(const vertex& v, T)
    S.erase(v);
  return (RELOAD_VBO);
}

// Balloon algorithm
uint ImgData::processBalloonSegmentation(HVecUS& data, Point12f& processParms)
{
  try {
    balloonModel.setStartRadius(processParms[2]);
    balloonModel.setMaxEdgeLen(processParms[3]);
    balloonModel.setStepSize(processParms[4]);
    balloonModel.setV0(processParms[5]);
    balloonModel.setI0(processParms[6]);
    balloonModel.setImageStd(processParms[7]);
    balloonModel.setVolumeStd(processParms[8]);
    balloonModel.setAreaStd(processParms[9]);
    balloonModel.setCurveStd(processParms[10]);

    int maxNbIter = (int)processParms[11];
    int maxNbRelax = maxNbIter / 4;

    balloonModel.init();

    Progress::start("Balloon algorithm", maxNbIter + maxNbRelax);

    QApplication::restoreOverrideCursor();
    int nIters = 0;
    while(balloonModel.step() > 0.005 && nIters < maxNbIter) {
      emit viewerUpdate();
      nIters++;
      if(!Progress::advance(nIters))
        return (RELOAD_VBO);
    }
    double volpen = balloonModel.getVolumePenalty();
    balloonModel.setVolumePenalty(0.0);
    nIters = 0;
    while(balloonModel.step() > 0.005 && nIters < maxNbRelax) {
      emit viewerUpdate();
      nIters++;
      if(!Progress::advance(maxNbIter + nIters))
        return (RELOAD_VBO);
    }
    balloonModel.setVolumePenalty(volpen);
  }
  catch(...) {
    SETSTATUS("Error during balloon segmentation.");
  }
  return (RELOAD_VBO);
}

// Convert mesh to cells
uint ImgData::processConvertCells(float cellMaxWall)
{
  // First check for bad corners
  forall(vertex v, S) {
    if(v->label == 0) {
      SETSTATUS("makeConvertCells::Error:Cannot make cells, unlabeled vertex found for " << Section);
      return (RELOAD_NONE);
    }
    forall(vertex n, S.neighbors(v)) {
      // get unique triangles
      vertex m = S.nextTo(v, n);
      if(!uniqueTri(S, v, n, m))
        continue;
      if(v->label < 0 and n->label < 0 and m->label < 0) {
        SETSTATUS("processConvertCells::Error:Cannot make cells, border triangles found for "
                  << Section << " (use fix corners)");
        return (RELOAD_NONE);
      }
    }
  }

  // Make centers
  typedef std::map<int, vertex> VMap;
  typedef pair<int, vertex> VMPair;
  VMap C;
  forall(const vertex& v, S) {
    forall(const vertex& n, S.neighbors(v)) {
      vertex m = S.nextTo(v, n);
      // get unique triangles
      if(!uniqueTri(S, v, n, m))
        continue;
      int label = getLabel(v, n, m);
      if(label <= 0)
        continue;
      // If new label add to map
      VMap::iterator i = C.find(label);
      if(i == C.end()) {
        vertex c;
        c->label = label;
        c->pos = 0;
        c->area = 0;
        i = C.insert(VMPair(label, c)).first;
      }
      vertex c = i->second;
      float area = triangleArea(v->pos, n->pos, m->pos);
      c->pos += (v->pos + n->pos + m->pos) / 3.0 * area;
      c->area += area;
    }
  }

  // Backup original graph and clear
  vvgraph T;
  T = S;
  S.clear();

  // Go through graph and count number of unqiue labels in neighborhood of each border vertex
  forall(const vertex& v, T) {
    v->type = ' ';
    if(v->label == -1) {
      std::set<int> L;
      forall(const vertex& n, T.neighbors(v))
        if(n->label > 0)
          L.insert(n->label);
      v->labcount = L.size();
    } else
      v->labcount = 0;
  }

  // Find all junctions
  for(VMap::iterator i = C.begin(); i != C.end(); i++) {
    vertex c = i->second;

    // Find a vertex labeled the same as c connected to an edge
    vertex v(0), n(0);
    bool found = false;
    forall(const vertex& w, T) {
      if(w->label == -1)
        forall(const vertex& m, T.neighbors(w))
          if(m->label == c->label) {
            found = true;
            v = w;
            n = m;
            break;
          }
      if(found)
        break;
    }
    if(!found)
      SETSTATUS("processConvertCells::Error:Label not found finding junctions");

    // Walk along edge until junction found
    vertex start = v;
    double length = 0;
    do {
      if(v->labcount > 2 or v->type == 'j' or length > worldToAbstract(cellMaxWall)) {       // at a junction
        length = 0;
        v->type = 'j';
      }
      n = T.prevTo(v, n);
      while(n->label != -1)
        n = T.prevTo(v, n);
      vertex m = n;
      length += (v->pos - n->pos).norm();
      n = v;
      v = m;
    } while(v != start);
  }

  // Connect centers to junctions (> 2 labels or margin and > 1 label)
  for(VMap::iterator i = C.begin(); i != C.end(); i++) {
    vertex c = i->second;

    // Find a vertex labeled the same as c connected to an edge
    vertex v(0), n(0);
    bool found = false;
    forall(const vertex& w, T) {
      if(w->label == -1)
        forall(const vertex& m, T.neighbors(w))
          if(m->label == c->label) {
            found = true;
            v = w;
            n = m;
            break;
          }
      if(found)
        break;
    }
    if(!found)
      SETSTATUS("processConvertCells::Error:Label not found connecting centers");

    // Walk along edge until junction found
    c->type = 'c';
    c->pos /= c->area;
    S.insert(c);
    vertex start = v;
    vertex pv(0);
    do {
      if(v->type == 'j') {       // at a junction
        v->type = 'j';
        if(!S.contains(v))
          S.insert(v);
        // Join centers to junctions
        if(S.valence(c) == 0)
          S.insertEdge(c, v);
        else
          S.spliceAfter(c, pv, v);
        pv = v;
      }
      n = T.prevTo(v, n);
      while(n->label != -1)
        n = T.prevTo(v, n);
      vertex m = n;
      n = v;
      v = m;
    } while(v != start);
  }

  // Join junctions to centers
  forall(const vertex& v, T) {
    if(v->type == 'j') {     // at a junction
      int labcount = v->labcount;
      vertex n = T.anyIn(v);
      vertex pn(0);
      while(labcount--) {
        while(n->label == -1)
          n = T.nextTo(v, n);
        int label = n->label;
        VMap::iterator i = C.find(label);
        if(i == C.end()) {
          SETSTATUS("processConvertCells::Error:Label not found connecting junctions, label:" << label);
          continue;
        }
        if(S.valence(v) == 0)
          S.insertEdge(v, i->second);
        else
          S.spliceAfter(v, pn, i->second);
        pn = i->second;
        while(n->label == label)
          n = T.nextTo(v, n);
      }
    }
  }

  // Remove points that are too close to each other
  std::set<vertex> D;
  forall(const vertex& c, S) {
    if(c->type != 'c')
      continue;
    forall(const vertex& n, S.neighbors(c)) {
      vertex m = S.nextTo(c, n);
      if((m->pos - n->pos).norm() < worldToAbstract(cellMaxWall) / 2.0) {
        if(D.find(n) != D.end() or D.find(m) != D.end())
          continue;
        if((n->margin and n->labcount == 1)or (!n->margin and n->labcount == 2))
          D.insert(n);
        else if((m->margin and m->labcount == 1)or (!m->margin and m->labcount == 2))
          D.insert(m);
      }
    }
  }
  forall(const vertex& v, D)
    if(S.contains(v))
      S.erase(v);

  // Connect junctions to each other
  forall(const vertex& c, S) {
    if(c->type != 'c')
      continue;
    forall(const vertex& n, S.neighbors(c)) {
      vertex m = S.nextTo(c, n);
      if(!S.edge(n, m))
        S.spliceBefore(n, c, m);
      if(!S.edge(m, n))
        S.spliceAfter(m, c, n);
    }
  }

  cells = true;
  // NextLabel = 0;
  SETSTATUS(Section << " make cells, " << S.size() << " vertices in mesh");
  return (RELOAD_VBO);
}
