//
// 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.
// 
#include <QFileDialog>

#include <StackProcess.hpp>
#include <Image.hpp>
#include <UnorderedMap.hpp>
#include <AutoTrim.hpp>
#include <GraphUtils.hpp> // nearestPointOnMesh


#include <algorithm>
#include <cfloat>

#include <CImg.h>

using namespace cimg_library;

typedef CImg<ushort> CImgUS;

#ifdef _MSC_VER
static float roundf(float v) // Defined to work only for positive values
{
  return floor(v + 0.5);
}
#endif

namespace mgx 
{
  // Get triangle list for cuda operations
  static void getTriangleList(Mesh *mesh, HVec3F& pts, HVec3U& tris, HVec3F& nrmls, int label)
  {
    // Cast because mesh attributes are maps (and lookup isn't const)
    const vvGraph &S = mesh->graph();
  
    const IntIntAttr &parents = mesh->parents();
  
    // Get triangle list
    uint tricount = 0;
    const std::vector<vertex> &av = mesh->activeVertices();
    bool has_selection = av.size() != S.size();
    std::set<vertex> vSet;
    forall(const vertex& v, av) {
      forall(const vertex& n, S.neighbors(v)) {
        vertex m = S.nextTo(v, n);
        if(has_selection and (!n->selected or !m->selected))
          continue;
        if(label >= 0 and label != mesh->getLabel(v, n, m, parents))
          continue;
        if(S.uniqueTri(v, n, m))
          tricount++;
  
        // Add vertices to set
        vSet.insert(v);
        vSet.insert(n);
        vSet.insert(m);
      }
    }
  
    // Put points in vector
    uint idx = 0;
    pts.resize(vSet.size());
    forall(const vertex& v, vSet) {
      pts[idx] = Point3f(v->pos);
      v->saveId = idx++;
    }
  
    idx = 0;
    tris.resize(tricount);
    nrmls.resize(tricount);
    forall(const vertex& v, av)
      forall(const vertex& n, S.neighbors(v)) {
        vertex m = S.nextTo(v, n);
        if(has_selection and (!n->selected or !m->selected))
          continue;
        if(label >= 0 and label != mesh->getLabel(v, n, m, parents))
          continue;
        if(!S.uniqueTri(v, n, m))
          continue;
        tris[idx] = Point3u(v->saveId, n->saveId, m->saveId);
        nrmls[idx++] = Point3f(((n->pos - v->pos) ^ (m->pos - v->pos)).normalize());
      }
  }
  
  // Get triangle list for cuda operations
  static void getPointList(const Mesh* mesh, HVec3F& pts, HVec3F& nrmls)
  {
    // Get triangle list
    uint idx = 0;
    const std::vector<vertex>& av = mesh->activeVertices();
    pts.resize(av.size());
    nrmls.resize(av.size());
    forall(const vertex& v, av) {
      pts[idx] = Point3f(v->pos);
      nrmls[idx] = Point3f(v->nrml);
      idx++;
    }
  }
  
  bool Annihilate::run(const Store* input, Store* output, const Mesh* mesh, float minDist, float maxDist,
                              bool fill, uint fillval)
  {
    HVecUS tdata(input->data());
  
    const Stack* stack = output->stack();
  
    // Data for cuda
    HVec3F pts;     // Points for triangles
    HVec3F nrmls;   // Normals for triangles
  
    // Get point list, nrmls, and bounding box
    getPointList(mesh, pts, nrmls);
    BoundingBox3f bBoxW;   // Bounding box in world coordinates
    forall(const Point3f& p, pts)
      bBoxW |= p;
  
    float distance = std::max(0.0f, std::max(-minDist, -maxDist));   // negative distance is outside mesh, expand bBox
    if(distance > 0) {
      Point3f distance3f(distance, distance, distance);
      bBoxW[0] -= distance3f;
      bBoxW[1] += distance3f;
    }
    BoundingBox3i bBoxI(stack->worldToImagei(bBoxW[0]), stack->worldToImagei(bBoxW[1]));
    Point3i imgSize(stack->size());
    bBoxI &= BoundingBox3i(Point3i(0, 0, 0), imgSize);
    Point3i size(bBoxI.size());
  
    // Buffer for mask
    HVecUS mask(size_t(size.x()) * size.y() * size.z(), 0);
  
    progressStart(QString("Annihilating Stack %1").arg(stack->userId()), 100);
  
    // call cuda to get stack mask
    if(nearMeshGPU(bBoxI[0], size, stack->step(), stack->origin(), pts, nrmls, minDist, maxDist, mask))
      return false;
    if(!progressAdvance(0))
      userCancel();
  
    size_t xIdx = 0;
    for(int z = 0; z < imgSize.z(); ++z)
      for(int y = 0; y < imgSize.y(); ++y)
        for(int x = 0; x < imgSize.x(); ++x, ++xIdx) {
          // Inside bBox
          if(x >= bBoxI[0].x() and x < bBoxI[1].x() and y >= bBoxI[0].y() and y < bBoxI[1].y()
             and z >= bBoxI[0].z() and z < bBoxI[1].z()) {
            size_t mIdx = ((size_t(z) - bBoxI[0].z()) * size.y() + y - bBoxI[0].y()) * size.x() + x
              - bBoxI[0].x();
            if(mask[mIdx]) {
              if(fill)
                output->data()[xIdx] = fillval;
              else
                output->data()[xIdx] = input->data()[xIdx];
              continue;
            }
          }
          if(fill and fillval == 0)
            output->data()[xIdx] = input->data()[xIdx];
          else
            output->data()[xIdx] = 0;
        }
    output->copyMetaData(input);
    output->changed();
    return true;
  }
  
  REGISTER_PROCESS(Annihilate);
  
  bool AutoScaleStack::run(const Store* input, Store* output)
  {
    unsigned int low = 65536, high = 0;
    const HVecUS& idata = input->data();
  
    for(size_t i = 0; i < idata.size(); ++i) {
      if(idata[i] < low)
        low = idata[i];
      if(idata[i] > high)
        high = idata[i];
    }
  
    unsigned int delta = high - low;
    if(delta == 0)
      delta = 65535u;
    HVecUS& data = output->data();
    if(delta == 65535u) {
      if(output != input) {
        data = idata;       // Nothing to scale
        output->changed();
      }
      return true;
    } else if(output != input) {
      for(size_t i = 0; i < idata.size(); ++i) {
        data[i] = (((unsigned int)(idata[i]) - low) * 65535u) / delta;
      }
    } else {
      for(size_t i = 0; i < data.size(); ++i) {
        data[i] = (((unsigned int)(data[i]) - low) * 65535u) / delta;
      }
    }
    output->copyMetaData(input);
    output->changed();
    return true;
  }
  
  REGISTER_PROCESS(AutoScaleStack);
  
  bool ApplyTransferFunction::run(Store* input, Store* output, float red, float green, float blue, float alpha)
  {
    TransferFunction::Colorf minColor, maxColor;
    const TransferFunction& fct = input->transferFct();
    float vmin = HUGE_VAL, vmax = -HUGE_VAL;
    if(input != output)
      output->allocate();
    HVecUS& data = input->data();
    HVecUS& outputData = output->data();
    float sum = red + green + blue + alpha;
    red /= sum;
    green /= sum;
    blue /= sum;
    alpha /= sum;
    for(size_t i = 0; i < data.size(); i++) {
      float f = float(data[i]) / 65535.0f;
      TransferFunction::Colorf c = fct.rgba(f);
      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(roundf(value * 65535));
      outputData[i] = (val > 65535 ? 65535 : val);
    }
    TransferFunction new_fct;
    new_fct.add_rgba_point(0, minColor);
    new_fct.add_rgba_point(1, maxColor);
    output->setTransferFct(new_fct);
    output->changed();
    output->copyMetaData(input);
    return true;
  }
  
  REGISTER_PROCESS(ApplyTransferFunction);
  
  typedef std::pair<Point3u, float> PixPair;
  static bool pixCompare(const PixPair& i, const PixPair& j) {
    return (i.second > j.second);
  }
  bool BlobDetect::run(const Store* input, Store* output, bool watershed, uint startlabel)
  {
    uint label = startlabel;
    HVecUS bdata = input->data();
    int id = input->stack()->userId();
    const Stack* stack = input->stack();
  
    // Sort pixels
    progressStart(QString("Stack %1 Blob Detect, sorting pixels").arg(id), 3);
    if(!progressAdvance(1))
      userCancel();
  
    Point3u size = stack->size();
    size_t SizeXYZ = size_t(size.x()) * size.y() * size.z();
    std::vector<PixPair> Pix(SizeXYZ);
    uint i = 0;
    for(uint z = 0; z < size.z(); z++)
      for(uint y = 0; y < size.y(); y++)
        for(uint x = 0; x < size.x(); x++) {
          Pix[i] = PixPair(Point3u(x, y, z), bdata[i]);
          i++;
        }
    if(!progressAdvance(2))
      userCancel();
    std::sort(Pix.begin(), Pix.end(), pixCompare);
  
    // Prepare result area
    output->reset();
    HVecUS& workData = output->data();
  
    progressStart(QString("Stack %1 Blob Detect, finding regions").arg(id), bdata.size());
    for(size_t i = 0; i < bdata.size(); i++) {
      // Find pixel offset
      const Point3u& p = Pix[i].first;
      size_t off = stack->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 >= (int)size.x())
        xe = (int)size.x() - 1;
  
      int ys = int(p.y() - 1);
      int ye = int(p.y() + 1);
      if(ys < 0)
        ys = 0;
      if(ye >= (int)size.y())
        ye = (int)size.y() - 1;
  
      int zs = int(p.z() - 1);
      int ze = int(p.z() + 1);
      if(zs < 0)
        zs = 0;
      if(ze >= (int)size.z())
        ze = (int)size.z() - 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 == (int)p.x() and y == (int)p.y() and z == (int)p.z())
              continue;
            size_t nboff = stack->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(QString("Blob Detect:too many labels"));
      if(i % 100000 == 0) {
        if(!progressAdvance(i))
          userCancel();
        updateState();
        updateViewer();
      }
    }
    // Clear all the background pixels
    for(size_t i = 0; i < bdata.size(); i++)
      if(workData[i] == 0xFFFF)
        workData[i] = 0;
  
    output->copyMetaData(input);
    output->changed();
  
    SETSTATUS("Stack " << id << " " << label - startlabel << " blobs found");
    return true;
  }
  
  REGISTER_PROCESS(BlobDetect);
  
  bool ClearWorkStack::run(Stack* stack, uint fillval)
  {
    Store* work = stack->work();
    work->allocate();
    if(fillval > 0xFFFF)
      fillval = 0xFFFF;
    HVecUS& data = work->data();
    std::fill_n(data.begin(), data.size(), ushort(fillval));
    work->changed();
    work->setFile("");
    return true;
  }
  
  REGISTER_PROCESS(ClearWorkStack);
  
  bool ClearMainStack::run(Stack* stack, uint fillval)
  {
    Store* main = stack->main();
    main->allocate();
    if(fillval > 0xFFFF)
      fillval = 0xFFFF;
    HVecUS& data = main->data();
    std::fill_n(data.begin(), data.size(), ushort(fillval));
    main->changed();
    main->setFile("");
    return true;
  }
  
  REGISTER_PROCESS(ClearMainStack);
  
  bool ClipStack::run(const Store* input, Store* output, bool trimStack)
  {
    // Setup matrices and planes for clip test
    Point3u clipDo;
    const Stack* stack = output->stack();
  
    // Indicator if in use
    clipDo[0] = (clip1()->enabled() or clip1()->grid() ? 1 : 0);
    clipDo[1] = (clip2()->enabled() or clip2()->grid() ? 1 : 0);
    clipDo[2] = (clip3()->enabled() or clip3()->grid() ? 1 : 0);
  
    // Frame matrix
    Matrix4d frm(stack->getFrame().worldMatrix());
  
    // Clip plane matrices
    Matrix4d clm[3];
    clm[0] = Matrix4d(clip1()->frame().inverse().worldMatrix());
    clm[1] = Matrix4d(clip2()->frame().inverse().worldMatrix());
    clm[2] = Matrix4d(clip3()->frame().inverse().worldMatrix());
  
    // Clipping planes, point normal form
    Point4f pn[6];
    pn[0] = Point4f(frm * clm[0]
                    * Point4d(clip1()->normal().x(), clip1()->normal().y(), clip1()->normal().z(), clip1()->width()));
    pn[1] = Point4f(frm * clm[0] 
                    * Point4d(-clip1()->normal().x(), -clip1()->normal().y(), -clip1()->normal().z(), clip1()->width()));
    pn[2] = Point4f(frm * clm[1]
                    * Point4d(clip2()->normal().x(), clip2()->normal().y(), clip2()->normal().z(), clip2()->width()));
    pn[3] = Point4f(frm * clm[1] 
                    * Point4d(-clip2()->normal().x(), -clip2()->normal().y(), -clip2()->normal().z(), clip2()->width()));
    pn[4] = Point4f(frm * clm[2]
                    * Point4d(clip3()->normal().x(), clip3()->normal().y(), clip3()->normal().z(), clip3()->width()));
    pn[5] = Point4f(frm * clm[2] 
                    * Point4d(-clip3()->normal().x(), -clip3()->normal().y(), -clip3()->normal().z(), clip3()->width()));
  
    // Put in host vector
    HVec4F Hpn;
    Hpn.resize(6);
    for(int i = 0; i < 6; i++)
      Hpn[i] = pn[i];
  
    // Call cuda to clip
    clipStackGPU(Point3i(stack->size()), stack->step(), stack->origin(), clipDo, Hpn, input->data(), output->data());

		// Call autotrim
	  if(trimStack){
	  	AutoTrim trim(*this);
		  bool result = trim.run(currentStack(), currentStack()->currentStore(), 0);
		  if(!result)
	  	  return false;
		}
  
    SETSTATUS("Stack " << stack->userId() << " clipped to intersection of clipping planes");
    output->copyMetaData(input);
    output->changed();
    return true;
  }
  
  REGISTER_PROCESS(ClipStack);
  
  bool CopyMainToWork::run(Stack* stack)
  {
    stack->work()->data() = stack->main()->data();
    stack->work()->changed();
    stack->work()->copyMetaData(stack->main());
    return true;
  }
  
  REGISTER_PROCESS(CopyMainToWork);
  
  bool CopyWorkToMain::run(Stack* stack)
  {
    stack->main()->data() = stack->work()->data();
    stack->main()->changed();
    stack->main()->copyMetaData(stack->work());
    return true;
  }
  
  REGISTER_PROCESS(CopyWorkToMain);
  
  bool StackMeshProcess::run(const Store* input, Store* output, Mesh* mesh,
                                                  int label, bool fill, uint fillval)
  {
    const Stack* stack = input->stack();
  
    if(fill)
      SETSTATUS("Filling stack " << stack->userId() << " from mesh");
    else
      SETSTATUS("Trimming stack " << stack->userId() << " from mesh");
  
    // Data for cuda
    HVec3F pts;     // Points for triangles
    HVec3U tris;    // Triangle vertices
    HVec3F nrmls;   // Normals for triangles
  
    // Get triangle list and bounding box
    getTriangleList(mesh, pts, tris, nrmls, label);
    if(pts.size() == 0)
      return false;
    BoundingBox3f bBoxW;
    forall(const Point3f& p, pts)
      bBoxW |= p;
  
    BoundingBox3i bBoxI(stack->worldToImagei(bBoxW[0]), 
                              stack->worldToImagei(bBoxW[1]) + Point3i(1, 1, 1));
    Point3i imgSize(stack->size());
    bBoxI &= BoundingBox3i(Point3i(0, 0, 0), imgSize);
    Point3i size = bBoxI.size();
  
    // Buffer for mask
    HVecUS mask(size_t(size.x()) * size.y() * size.z(), 0);
  
    // call cuda to get stack mask
    progressStart(QString("Trimming Stack %1").arg(stack->userId()), 100);
    if(insideMeshGPU(bBoxI[0], size, stack->step(), stack->origin(), pts, tris, nrmls, mask))
      return false;
    if(!progressAdvance(0))
      userCancel();
  
    size_t xIdx = 0;
    for(int z = 0; z < imgSize.z(); ++z)
      for(int y = 0; y < imgSize.y(); ++y)
        for(int x = 0; x < imgSize.x(); ++x, ++xIdx) {
          // Inside bBox
          if(x >= bBoxI[0].x() and x < bBoxI[1].x() and y >= bBoxI[0].y() and y < bBoxI[1].y()
             and z >= bBoxI[0].z() and z < bBoxI[1].z()) {
            size_t mIdx = ((size_t(z) - bBoxI[0].z()) * size.y() + y - bBoxI[0].y()) * size.x() + x
              - bBoxI[0].x();
            if(mask[mIdx]) {
              if(fill)
                output->data()[xIdx] = fillval;
              else
                output->data()[xIdx] = input->data()[xIdx];
              continue;
            }
          }
          if(fill and fillval == 0)
            output->data()[xIdx] = input->data()[xIdx];
          else
            output->data()[xIdx] = 0;
        }
  
    output->copyMetaData(input);
    output->changed();
    return true;
  }
  
  REGISTER_PROCESS(FillStackToMesh);
  REGISTER_PROCESS(TrimStackProcess);
  
  bool FillStack3D::run(const Store* input, Store* output, Mesh* mesh)
  {
    const Stack* stack = input->stack();
  
    SETSTATUS("Filling stack " << stack->userId() << " from 3D mesh");
  
    // Data for cuda
    HVec3F pts;     // Points for triangles
    HVec3U tris;    // Triangle vertices
    HVec3F nrmls;   // Normals for triangles
                    //    Point3f bBoxW[2];  // Bounding box in world coordinates
  
    // Get list of labels
    std::set<int> Labels;
    forall(const vertex& v, mesh->graph())
      if(v->label > 0)
        Labels.insert(v->label);
  
    // Loop over all labels
    progressStart(QString("Filling Stack %1 from 3D mesh").arg(stack->userId()), 0);
    forall(int label, Labels) {
      std::cout << "label: " << label << std::endl;
      // Get triangle list and bounding box
      getTriangleList(mesh, pts, tris, nrmls, label);
      BoundingBox3f bBoxW;
      forall(const Point3f& p, pts)
        bBoxW |= p;
      BoundingBox3i bBoxI(stack->worldToImagei(bBoxW[0]), stack->worldToImagei(bBoxW[1]) + Point3i(1, 1, 1));
      Point3i imgSize(stack->size());
      bBoxI &= BoundingBox3i(Point3i(0, 0, 0), imgSize);
      Point3i size = bBoxI.size();
      // Buffer for mask
      HVecUS mask(size_t(size.x()) * size.y() * size.z(), 0);
  
      // call cuda to get stack mask
      if(insideMeshGPU(bBoxI[0], size, stack->step(), stack->origin(), pts, tris, nrmls, mask))
        return false;
      if(!progressAdvance(0))
        userCancel();
      for(int z = bBoxI[0].z(); z < bBoxI[1].z(); ++z)
        for(int y = bBoxI[0].y(); y < bBoxI[1].y(); ++y) {
          size_t xIdx = stack->offset(bBoxI[0].x(), y, z);
          for(int x = bBoxI[0].x(); x < bBoxI[1].x(); ++x, ++xIdx) {
            size_t mIdx = ((size_t(z) - bBoxI[0].z()) * size.y() + y - bBoxI[0].y()) * size.x() + x
              - bBoxI[0].x();
            if(mask[mIdx])
              output->data()[xIdx] = label;
          }
        }
    }
    output->copyMetaData(input);
    output->changed();
    return true;
  }
  REGISTER_PROCESS(FillStack3D);
  
  bool SwapStacks::run(Stack* stack)
  {
    if(!stack->main()) {
      setErrorMessage("Error, the main stack is not initialized");
      return false;
    }
    if(!stack->work()) {
      setErrorMessage("Error, the work stack is not initialized");
      return false;
    }
    using std::swap;
    swap(stack->work()->data(), stack->main()->data());
    swapMetaData(stack->work(), stack->main());
  
    stack->main()->changed();
    stack->work()->changed();
    return true;
  }
  
  REGISTER_PROCESS(SwapStacks);
  
  bool CopySwapStacks::run(const QString& storeStr, const QString& actionStr)
  {
    Stack* stack1 = stack(0);
    if(!stack1) {
      setErrorMessage("Error, bad action, stack 1 empty");
      return false;
    }
    Stack* stack2 = stack(1);
    if(!stack2) {
      setErrorMessage("Error, bad action, stack 2 empty");
      return false;
    }
    Store* store1 = (storeStr == "Main" ? stack1->main() : stack1->work());
    if(!store1) {
      setErrorMessage("Error, bad action, store 1 empty");
      return false;
    }
    Store* store2 = (storeStr == "Main" ? stack2->main() : stack2->work());
    if(!store2) {
      setErrorMessage("Error, bad action, store 2 empty");
      return false;
    }
  
    bool diffStep = (stack1->step() != stack2->step());
    bool diffOrigin = (stack1->origin() != stack2->origin());
    bool diffSize = (stack1->size() != stack2->size());
  
    if(actionStr == "1 -> 2") {
      if(stack2->storeSize() == 0) {
        stack2->setStep(stack1->step());
        stack2->setOrigin(stack1->origin());
        stack2->setSize(stack1->size());
      } else if(diffStep or diffOrigin or diffSize)
        return setErrorMessage("Error, Stack 1 and 2 have different size, step, or origin.");
      store2->data() = store1->data();
      store2->copyMetaData(store1);
      store2->changed();
      return true;
    } else if(actionStr == "1 <- 2") {
      if(stack1->storeSize() == 0) {
        stack1->setStep(stack2->step());
        stack1->setOrigin(stack2->origin());
        stack1->setSize(stack2->size());
      } else if(diffStep or diffOrigin or diffSize) {
        setErrorMessage("Error, Stack 1 and 2 have different size, step, or origin.");
        return false;
      }
      store1->data() = store2->data();
      store1->copyMetaData(store2);
      store1->changed();
      return true;
    } else if(actionStr == "1 <-> 2") {
      if(diffStep or diffOrigin or diffSize) {
        setErrorMessage("Error, Stack 1 and 2 have different size, step, or origin.");
        return false;
      }
      HVecUS tempData(stack1->main()->data());
      Point3f tempStep(stack1->step());
      Point3f tempOrigin(stack1->origin());
      Point3u tempSize(stack1->size());
  
      bool l1 = store1->labels();
      bool l2 = store2->labels();
      QString f1 = store1->file();
      QString f2 = store2->file();
  
      stack1->setStep(stack2->step());
      stack1->setOrigin(stack2->origin());
      stack1->setSize(stack2->size());
      store1->data() = store2->data();
  
      stack2->setStep(tempStep);
      stack2->setOrigin(tempOrigin);
      stack2->setSize(tempSize);
      store2->data() = tempData;
  
      store1->setFile(f2);
      store1->setLabels(l2);
      store2->setFile(f1);
      store2->setLabels(l1);
  
      store1->changed();
      store2->changed();
      return (true);
    } else {
      setErrorMessage("Error, bad action:" + actionStr);
      return false;
    }
    return true;
  }
  
  REGISTER_PROCESS(CopySwapStacks);
  
  bool ResizeCanvas::run(Stack* stack, bool isRelative, bool center, Point3i ds)
  {
    Point3i oldSize(stack->size());
    Point3i newSize(stack->size());
  
    HVecUS& MainData = stack->main()->data();
    HVecUS& WorkData = stack->work()->data();
  
    if(isRelative)
      newSize += ds;
    else
      for(int i = 0; i < 3; ++i)
        if(ds[i] > 0)
          newSize[i] = ds[i];

    if(newSize == oldSize)
      return true;

    // Check new size
    for(int i = 0; i < 3; ++i)
      if(newSize[i] < 0)
        throw(QString("Negative size"));
  
    {
      HVecUS mi = resize(MainData, oldSize, newSize, center);
      MainData.swap(mi);
    }
    {
      HVecUS wi = resize(WorkData, oldSize, newSize, center);
      WorkData.swap(wi);
    }
    stack->setSize(Point3u(newSize));
  
    stack->main()->changed();
    stack->work()->changed();
  
    return true;
  }
  
  REGISTER_PROCESS(ResizeCanvas);
  
  bool ScaleStack::run(Stack* stack, Point3f newsize, bool percent)
  {
    if(newsize.x() <= 0 and newsize.y() <= 0.0 and newsize.z() <= 0) {
      setErrorMessage("Error, no rescale specified.");
      return false;
    }
    Point3u size = stack->size();
    Point3f step = stack->step();
  
    HVecUS& MainData = stack->main()->data();
    HVecUS& WorkData = stack->work()->data();
  
    CImgUS mainImage(MainData.data(), size.x(), size.y(), size.z(), 1, false);
    CImgUS workImage(WorkData.data(), size.x(), size.y(), size.z(), 1, false);
  
    if(percent)
      for(int i = 0; i < 3; i++)
        newsize[i] *= size[i] / 100.0;
  
    for(int i = 0; i < 3; i++)
      if(newsize[i] > 0) {
        step[i] *= size[i] / newsize[i];
        size[i] = newsize[i];
      }
  
    mainImage.resize(size.x(), size.y(), size.z(), 1, (stack->main()->labels() ? 1 : 5));
    workImage.resize(size.x(), size.y(), size.z(), 1, (stack->work()->labels() ? 1 : 5));
  
    stack->setSize(size);
    stack->setStep(step);
  
    memcpy(MainData.data(), mainImage.data(), MainData.size() * 2);
    memcpy(WorkData.data(), workImage.data(), WorkData.size() * 2);
  
    stack->main()->changed();
    stack->work()->changed();
  
    return true;
  }
  
  REGISTER_PROCESS(ScaleStack);
  
  bool ShiftStack::run(Stack* stack, bool origin, Point3i ds)
  {
    Point3u size = stack->size();
  
    Store* main = stack->main();
    Store* work = stack->work();
  
    HVecUS& MainData = main->data();
    HVecUS& WorkData = work->data();
  
    if(origin)
      stack->setOrigin(stack->origin() + Point3f(ds)); // RSS why is ds Point3i?
    else {
      CImgUS mainImage(MainData.data(), size.x(), size.y(), size.z(), 1, true);
      mainImage.shift(ds.x(), ds.y(), ds.z(), 0, 0);
  
      CImgUS workImage(WorkData.data(), size.x(), size.y(), size.z(), 1, true);
      workImage.shift(ds.x(), ds.y(), ds.z(), 0, 0);
    }
  
    main->changed();
    work->changed();
    return true;
  }
  
  REGISTER_PROCESS(ShiftStack);
  
  bool StackRelabel::run(Stack* stack, const Store* store, Store* output, int start, int step)
  {
    if(step < 1)
      step = 1;
    if(start < 1)
      start = 1;
    progressStart(QString("Relabeling current stack"), store->size());
    std::unordered_map<int, int> relabel_map;
    int prog = store->size() / 100;
    const HVecUS& data = store->data();
    HVecUS& odata = output->data();
    ushort label = 1;
    for(size_t i = 0; i < data.size(); ++i) {
      if(data[i] > 0u) {
        std::unordered_map<int, int>::iterator found = relabel_map.find(data[i]);
        if(found != relabel_map.end())
          odata[i] = found->second;
        else {
          relabel_map[data[i]] = label;
          odata[i] = label;
          label++;
        }
      }
      if((i % prog) == 0 and !progressAdvance(i))
        userCancel();
    }
    // Find a random arangement of labels
    std::vector<ushort> labels(label);
    for(ushort i = 0; i < (ushort)labels.size(); ++i)
      labels[i] = ushort(start + i * step);
    ushort last = labels.back();
    std::random_shuffle(labels.begin(), labels.end());
  #pragma omp parallel for
    for(size_t i = 0; i < data.size(); ++i) {
      if(odata[i] > 0u) {
        odata[i] = labels[odata[i] - 1];
      }
    }
  
    if(!progressAdvance(data.size()))
      userCancel();
    output->changed();
    output->copyMetaData(store);
    stack->setLabel(last);
    Information::out << "Next label = " << stack->viewLabel() << endl;
    SETSTATUS(QString("Total number of labels created = %1. Last label = %2").arg(relabel_map.size()).arg(last));
    return true;
  }
  
  REGISTER_PROCESS(StackRelabel);
  
  bool ChangeVoxelSize::run(Stack* stack, Point3f nv)
  {
    if(nv.x() <= 0)
      nv.x() = stack->step().x();
    if(nv.y() <= 0)
      nv.y() = stack->step().y();
    if(nv.z() <= 0)
      nv.z() = stack->step().z();
    Point3f ratio = nv / stack->step();
    stack->setStep(nv);
    stack->setOrigin(multiply(stack->origin(), ratio));
    stack->main()->changed();
    stack->work()->changed();
    return true;
  }
  
  REGISTER_PROCESS(ChangeVoxelSize);
  
  bool ReverseStack::run(Store* output, const Store* input, bool rx, bool ry, bool rz)
  {
    const Stack* stack = output->stack();
    Point3u size = stack->size();
    if(input == output) {
      HVecUS& data = output->data();
      if(rx and ry and rz) {
        ushort* s = &data[0];
        ushort* e = s + data.size() - 1;
        uint nb_iter = data.size() >> 1;
        for(uint i = 0; i < nb_iter; ++i) {
          std::swap(*e--, *s++);
        }
        output->changed();
      } else if(rx xor ry xor rz)     // Only one if true
      {
        int dx = 1, dy = 1, dz = 1;
        Point3u start;
        if(rx) {
          dx = -1;
          start.x() = size.x() - 1;
          size.x() >>= 1;
        } else if(ry) {
          dy = -1;
          start.y() = size.y() - 1;
          size.y() >>= 1;
        } else {
          dz = -1;
          start.z() = size.z() - 1;
          size.z() >>= 1;
        }
        for(uint z = 0, z1 = start.z(); z < size.z(); ++z, z1 += dz)
          for(uint y = 0, y1 = start.y(); y < size.y(); ++y, y1 += dy)
            for(uint x = 0, x1 = start.x(); x < size.x(); ++x, x1 += dx) {
              std::swap(data[stack->offset(x, y, z)], data[stack->offset(x1, y1, z1)]);
            }
        output->changed();
      } else if(rx and ry) {
        uint size_xy = size.x() * size.y();
        uint to_shift = size_xy >> 1;
        for(uint z = 0; z < size.z(); ++z) {
          ushort* s = &data[stack->offset(0, 0, z)];
          ushort* e = s + size_xy - 1;
          for(uint i = 0; i < to_shift; ++i) {
            std::swap(*e--, *s++);
          }
        }
        output->changed();
      } else if(rz)     // and rx or ry
      {
        Point3u start;
        int dx = 1, dy = 1;
        start.z() = size.z() - 1;
        if(rx) {
          start.x() = size.x() - 1;
          dx = -1;
        } else       // if(ry)
        {
          start.y() = size.y() - 1;
          dy = -1;
        }
        uint half_z = size.z() >> 1;
        for(uint z = 0, z1 = start.z(); z < half_z; ++z, --z1)
          for(uint y = 0, y1 = start.y(); y < size.y(); ++y, y1 += dy)
            for(uint x = 0, x1 = start.x(); x < size.x(); ++x, x1 += dx) {
              std::swap(data[stack->offset(x, y, z)], data[stack->offset(x1, y1, z1)]);
            }
        if(size.z() & 1)       // If there is an off number of planes
        {
          uint z = half_z;
          uint end_y = (ry ? size.y() >> 1 : size.y());
          uint end_x = (rx ? size.x() >> 1 : size.x());
          for(uint y = 0, y1 = start.y(); y < end_y; ++y, y1 += dy)
            for(uint x = 0, x1 = start.x(); x < end_x; ++x, x1 += dx) {
              std::swap(data[stack->offset(x, y, z)], data[stack->offset(x1, y1, z)]);
            }
        }
        output->changed();
      }
    } else if(rx or ry or rz) {
      const HVecUS& data = input->data();
      HVecUS& out = output->data();
      Point3u size = stack->size();
      int dx = (rx ? -1 : 1), dy = (ry ? -1 : 1), dz = (rz ? -1 : 1);
      uint startx = (rx ? size.x() - 1 : 0);
      uint starty = (ry ? size.y() - 1 : 0);
      uint startz = (rz ? size.z() - 1 : 0);
      for(uint z1 = 0, z2 = startz; z1 < size.z(); ++z1, z2 += dz)
        for(uint y1 = 0, y2 = starty; y1 < size.y(); ++y1, y2 += dy)
          for(uint x1 = 0, x2 = startx; x1 < size.x(); ++x1, x2 += dx)
            out[stack->offset(x2, y2, z2)] = data[stack->offset(x1, y1, z1)];
      output->copyMetaData(input);
      output->changed();
    }
    return true;
  }
  
  REGISTER_PROCESS(ReverseStack);
  
  bool StackRelabelFromMesh::run(Stack* stack, const Store* store, Store* output, const Mesh* mesh,
                                        bool delete_unknown)
  {
    VVGraphVec regions;
    VtxIntMap cells;
    int n = mesh->getConnectedRegions(mesh->graph(), regions, cells);
    SETSTATUS("Number of cells = " << n);
    std::vector<Point3f> centers(regions.size());
    std::vector<ushort> labels(regions.size());
    std::map<ushort, ushort> label_mapping;
    Point3u stkSize = stack->size();
    const HVecUS& in_data = store->data();
    HVecUS& data = output->data();
    ushort max_label = 0;
  
    for(size_t i = 0; i < regions.size(); ++i) {
      const vvGraph& S = regions[i];
      short lab = ushort(S.any()->label & 0xFFFF);
      if(lab > max_label)
        max_label = lab;
      labels[i] = lab;
      Point3f& c = centers[i];
      forall(const vertex& v, S)
        c += Point3f(v->pos);
      c /= S.size();
      Point3i imgPos = stack->worldToImagei(c);
      if(imgPos.x() >= 0 and imgPos.y() >= 0 and imgPos.z() >= 0 and (size_t) imgPos.x() < stkSize.x()
         and (size_t) imgPos.y() < stkSize.y() and (size_t) imgPos.z() < stkSize.z()) {
        ushort l = in_data[stack->offset(imgPos)];
        label_mapping[l] = lab;
        SETSTATUS("Mapping " << l << " to " << lab);
      } else
        SETSTATUS("No mapping for cell labeled " << lab);
    }
  
    for(size_t i = 0; i < data.size(); ++i) {
      std::map<ushort, ushort>::const_iterator found = label_mapping.find(in_data[i]);
      if(found != label_mapping.end())
        data[i] = found->second;
      else if(delete_unknown)
        data[i] = 0;
      else {
        label_mapping[in_data[i]] = ++max_label;
        data[i] = max_label;
      }
    }
  
    output->copyMetaData(store);
    output->changed();
  
    return true;
  }
  REGISTER_PROCESS(StackRelabelFromMesh);

  bool DeleteOutsideMesh::run(Stack* stack, const Store* input, Store* output, const Mesh* mesh, double threshold)
  {
    std::vector<Point3d> centersMesh(65536, Point3d());
    std::vector<double> weightsMesh(65536, 0);
  
    // Calculate mesh centers
    const vvGraph &S = mesh->graph();
    forall(const vertex& v, S)
      forall(const vertex& n, S.neighbors(v)) {
        vertex m = S.nextTo(v, n);
        if(!S.uniqueTri(v, n, m))
          continue;
        // Calculate volume enclosed by mesh and save triangles
        int label = v->label;
        if(label <= 0)
          throw QString("%1::run Unlabeled vertex").arg(name());
        if(label != n->label or label != m->label)
          throw QString("%1::run Multiply labeled cell").arg(name());

        double volume = signedTetraVolume(v->pos, n->pos, m->pos);
        centersMesh[label] += volume * Point3d(v->pos + n->pos + m->pos) / 4.0;
        weightsMesh[label] += volume; 
      }

    #pragma omp parallel for
    for(size_t label = 1; label < 65536; label++)
      if(weightsMesh[label] <= 0)
        continue;
      else
        centersMesh[label] /= weightsMesh[label];

    const HVecUS& inputData = input->data();
    std::vector<Point3d> centersStack(65536, Point3d());
    std::vector<size_t> weightsStack(65536, 0);
    auto size = stack->size();
    for(uint z = 0; z < size.z(); z++)
      for(uint y = 0; y < size.y(); y++)
        for(uint x = 0; x < size.x(); x++) {
          int label = inputData[stack->offset(x,y,z)];
          if(label <= 0)
            continue;
          centersStack[label] += Point3d(stack->imageToWorld(Point3i(x, y, z)));
          weightsStack[label]++;
        }

    #pragma omp parallel for
    for(int label = 0; label < 65536; label++)
      if(weightsStack[label] > 0)
        centersStack[label] /= weightsStack[label];

    HVecUS& outputData = output->data();
    size_t totalSize = size_t(size.x()) * size.y() * size.z();
    #pragma omp parallel for
    for(size_t i = 0; i < totalSize; i++) {
      int label = inputData[i];
      if(label > 0 and weightsMesh[label] > 0 and weightsStack[label] > 0 and norm(centersMesh[label] - centersStack[label]) <= threshold)
        outputData[i] = label;
      else
        outputData[i] = 0;
    }

    output->copyMetaData(input);
    output->changed();
  
    return true;
  }
  REGISTER_PROCESS(DeleteOutsideMesh);
 
  bool SaveTransform::initialize(QWidget* parent)
  {
    QString filename = parm("Filename");
    if(filename.isEmpty())
      filename = QFileDialog::getSaveFileName(parent, "Choose transform file to save", 
                                          QDir::currentPath(), "Text files (*.txt)");
    if(filename.isEmpty())
      return false;
    if(!filename.endsWith(".txt", Qt::CaseInsensitive))
      filename += ".txt";
    setParm("Filename", filename);
    return true;
  }
  
  bool SaveTransform::run(Stack* stack, const QString& filename)
  {
    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly)) {
      setErrorMessage(QString("File '%1' cannot be opened for writing").arg(filename));
      return false;
    }
    QTextStream out(&file);
  
    // Write matrix to file
    Matrix4d m(stack->getFrame().worldMatrix());
    out << transpose(m) << endl;
    file.close();
  
    SETSTATUS(QString("Transform saved to file: %1").arg(filename));
    return true;
  }
  REGISTER_PROCESS(SaveTransform);
  
  bool LoadTransform::initialize(QWidget* parent)
  {
    QString filename = parm("Filename");
    if(filename.isEmpty())
      filename = QFileDialog::getOpenFileName(parent, "Choose transform file to load", QDir::currentPath(),
                                              "Text files (*.txt)");
    if(filename.isEmpty())
      return false;
    if(!filename.endsWith(".txt", Qt::CaseInsensitive))
      filename += ".txt";
    setParm("Filename", filename);
    return true;
  }
  
  bool LoadTransform::run(Stack* stack, const QString& filename)
  {
    QFile file(filename);
    if(!file.open(QIODevice::ReadOnly)) {
      setErrorMessage(QString("File '%1' cannot be opened for reading").arg(filename));
      return false;
    }
    QTextStream out(&file);
  
    // Read in transform
    Matrix4d m;
    QTextStream in(&file);
    in >> m;
    file.close();
  
    // Set frame or transform
    if(stack->showTrans())
      stack->trans().setFromMatrix((~m).c_data());
    else
      stack->frame().setFromMatrix((~m).c_data());
  
    SETSTATUS(QString("Transform loaded from file: %1").arg(filename));
    return true;
  }
  REGISTER_PROCESS(LoadTransform);


  bool DistanceStackToMesh::run(Stack* s1, Store* store1, Store* store2, Mesh* m, double mult)
  {
    HVecUS& data2 = store2->data();
    Point3i imgSize(s1->size());

    int xDim = imgSize.x();
    int yDim = imgSize.y();
    int zDim = imgSize.z();

    #pragma omp parallel for schedule(guided)
    for(int x=0; x<=xDim; x++){
      for(int y=0; y<=yDim; y++){
        for(int z=0; z<=zDim; z++){
          Point3i pOut (x,y,z);
          Point3d pOutDouble(s1->imageToWorld(pOut));
          //distance to mesh
          Point3d p = nearestPointOnMesh(pOutDouble, m->graph());
          double distance = norm(p - pOutDouble);
          distance*=mult;
          data2[s1->offset(pOut)] = (int)std::round(distance);
        }
      }
    }

    store2->changed();
    return true;
  }
  REGISTER_PROCESS(DistanceStackToMesh);

   std::set<int> findNeighborsLabel(Store* store, Point3i imgSize, Point3i p, int radius)
   {
     HVecUS& dataV = store->data();
     const Stack* stack = store->stack();
     std::set<int> result;
     Point3i pLocal = p;
     for(int z = -radius; z <= radius; ++z)
       for(int y = -radius; y <= radius; ++y)
         for(int x = -radius; x <= radius; ++x) {
           pLocal.x() = p.x() + x;
           pLocal.y() = p.y() + y;
           pLocal.z() = p.z() + z;
           if(pLocal.x()<0 or pLocal.x()>=imgSize.x() or pLocal.y()<0 or pLocal.y()>=imgSize.y() or pLocal.z()<0 or pLocal.z()>=imgSize.z()) continue;
           int label = dataV[stack->offset(pLocal)];
           result.insert(label);
         }

      return result;
   }



}
