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

#include <set>
#include <vector>

using std::cout;
using std::cerr;
using std::endl;
using std::string;

namespace mgx 
{
  // Color gradient along Z axis only
  struct colorGradOP 
  {
    const float div;
  
    Point3i size;
    Point3i stride;
    ushort* Src;
    ushort* Dst;
  
    colorGradOP(float _div, DeviceVus* _Src, DeviceVus* _Dst)
      : div(_div) , Src(devP(*_Src)) , Dst(devP(*_Dst)) {}
  
    void setSizes(const Point3i&, const Point3i& _size)
    {
      size = _size;
      stride = getStride(size);
    }
  
    __device__ __host__ void operator()(uint idx) const
    {
      int x, y;
      getXY(idx, size, x, y);
      float grad, nextgrad, prevgrad = 0;
  
      size_t zIdx = getOffset(x, y, 0, size);
      for(int z = 0; z < size.z(); ++z, zIdx += stride.z()) {
        if(z < size.z() - 1)
          nextgrad = float(Src[zIdx] - Src[zIdx + stride.z()]) / div;
        else
          nextgrad = 0;
        grad = (nextgrad + prevgrad) / div;
        prevgrad = nextgrad;
  
        grad *= grad;
        if(grad > 0xFFFF)
          grad = 0xFFFF;
        Dst[zIdx] = (ushort)(grad);
      }
    }
  };
  
  // Edge detect along Y axis only, image is flipped
  struct edgeDetectOP
  {
    const ushort low;
    const ushort high;
    const float mult;
    const ushort fillval;
  
    Point3i size;
    Point3i stride;
    ushort* Src;
    ushort* Dst;
  
    edgeDetectOP(ushort _low, ushort _high, float _mult, ushort _fillval, DeviceVus* _Src, DeviceVus* _Dst)
      : low(_low), high(_high), mult(_mult), fillval(_fillval), Src(devP(*_Src)), Dst(devP(*_Dst)) {}
  
    void setSizes(const Point3i&, const Point3i& _size)
    {
      size = _size;
      stride = getStride(size);
    }
  
    __device__ __host__ void operator()(uint idx) const
    {
      int x, z;
      getXZ(idx, size, x, z);
  
      int count = 0;
      float avg = 0;
      size_t yIdx = getOffset(x, 0, z, size);
      for(int y = 0; y < size.y(); ++y, yIdx += stride.y()) {
        uint s = Src[yIdx];
        if(s >= low) {
          count++;
          avg += s;
        }
      }
      int thresh = high;
      if(count > 0) {
        avg /= count;
        if(avg < high)
          thresh = avg;
      }
  
      int found = 0;
      yIdx = getOffset(x, size.y() - 1, z, size);
      for(int y = size.y() - 1; y >= 0; --y, yIdx -= stride.y()) {
        if(found == 0 && Src[yIdx] >= thresh)
          found = 1;
        if(thresh > 0 and found == 1)
          Dst[yIdx] = fillval;
        else
          Dst[yIdx] = 0;
      }
    }
  };
  
  // Clip stack to clip planes
  struct clipStackOP 
  {
    const Point3f step;
    const Point3f shift;
    const Point3u clipDo;
    Point4f* pn;
  
    Point3i base;
    Point3i size;
    Point3i stride;
    ushort* Src;
    ushort* Dst;
  
    clipStackOP(const Point3f& _step, const Point3f& _shift, const Point3u& _clipDo, DeviceV4f& _pn, DeviceVus* _Src, DeviceVus* _Dst)
      : step(_step), shift(_shift), clipDo(_clipDo), pn(devP(_pn)), Src(devP(*_Src)), Dst(devP(*_Dst)) {}
  
    void setSizes(const Point3i& _base, const Point3i& _size)
    {
      base = _base;
      size = _size;
      stride = getStride(size);
    }
  
    __device__ __host__ void operator()(uint idx) const
    {
      int x, y;
      getXY(idx, size, x, y);
  
      // Copy && clip
      size_t zIdx = getOffset(x, y, 0, size);
      for(int z = 0; z < size.z(); z++, zIdx += stride.z()) {
        Point3f wp = imageToWorld(base + Point3i(x, y, z), step, shift);
        if(clipPoint(wp, clipDo, &pn[0]))
          Dst[zIdx] = 0;
        else
          Dst[zIdx] = Src[zIdx];
      }
    }
  };
  
  // Calculate mask contained by a volume
  struct insideMeshOP 
  {
    const Point3f step;
    const Point3f shift;
    const Point3f* pts;
    const Point3u* tris;
    const Point3f* nrmls;
    const Point3i imgBase;
    const uint* gridIdx;
    const uint* gridCnt;
    const uint* gridTris;
    const Point3i gridSize;
    const Point3i gridPix;
  
    Point3i size;
    Point3i base;
    Point3i stride;
    ushort* Dst;
  
    insideMeshOP(const Point3f& _step, const Point3f& _shift, const Point3f* _pts, 
        const Point3u* _tris, const Point3f* _nrmls, const Point3i& _imgBase, const uint* _gridIdx, 
        const uint* _gridCnt, const uint* _gridTris, const Point3i _gridSize, 
        const Point3i& _gridPix, DeviceVus& _Dst)
      : step(_step), shift(_shift), pts(_pts), tris(_tris), nrmls(_nrmls), imgBase(_imgBase), 
        gridIdx(_gridIdx), gridCnt(_gridCnt), gridTris(_gridTris), gridSize(_gridSize), 
        gridPix(_gridPix), Dst(devP(_Dst)) {}
  
    void setSizes(const Point3i& _base, const Point3i& _size)
    {
      base = _base;
      size = _size;
      stride = getStride(size);
    }
  
    __device__ __host__ int rayXTriangleIntersect(const Point3f& pt, const Point3f& t0, const Point3f& t1,
                                                  const Point3f& t2, const Point3f& n, Point3f& intp)
    {
      // This routine is optimized for axis aligned rays in X direction
      if(n.x() == 0)     // triangle lies in yz plane
        return 0;
  
      Point3f w0 = pt - t0;
      float a = -n * w0;
      float r = a / n.x();
      if(r < 0.0f)      // ray goes away from triangle
        return 0;       // => no intersect
  
      intp = pt;
      intp.x() += r;     // intersect point of ray and plane
  
      // check if intersect point inside triangle
      Point3f u = t1 - t0;
      Point3f v = t2 - t0;
      float uu, uv, vv, wu, wv, d;
      uu = u * u;
      uv = u * v;
      vv = v * v;
      Point3f w = intp - t0;
      wu = w * u;
      wv = w * v;
      d = uv * uv - uu * vv;
  
      // get and test parametric coords
      float s, t;
      s = (uv * wv - vv * wu) / d;
      if(s < -FLT_EPSILON || s > 1.0f + FLT_EPSILON)     // outside
        return 0;
      t = (uv * wu - uu * wv) / d;
      if(t < -FLT_EPSILON || (s + t) > 1.0f + FLT_EPSILON)     // outside
        return 0;
  
      return 1;     // inside
    }
  
    __device__ __host__ void operator()(uint idx)
    {
      Point3f intp;
      int y, z;
      getYZ(idx, size, y, z);
  
      // Look for intersections down a column of pixels
      int x = 0;
      Point3f pt = imageToWorld(imgBase + base + Point3i(x, y, z), step, shift);
      pt.x() -= 100;
  
      // Grid information
      int yshift = base.y();
      int zshift = base.z();
      ;
      Point3i gridPos(0, (yshift + y) / gridPix.y(), (zshift + z) / gridPix.z());
      int idxG = gridPos.z() * gridSize.y() + gridPos.y();
      do {
        // Loop through triangle list to find closest triangle along x axis
        Point3f minpt;
        float mindist = FLT_MAX;
        bool inside = false;
        int idxT = gridIdx[idxG];
        for(uint i = 0; i < gridCnt[idxG]; ++i, ++idxT) {
          int idxTri = gridTris[idxT];
          Point3f t0 = pts[tris[idxTri].x()];
          Point3f t1 = pts[tris[idxTri].y()];
          Point3f t2 = pts[tris[idxTri].z()];
          Point3f n = nrmls[idxTri];
  
          // Skip if triangle is behind
          float maxx = max(t0.x(), max(t1.x(), t2.x()));
          if(maxx < pt.x())
            continue;
  
          // See if we already have a closer triangle
          float minx = min(t0.x(), min(t1.x(), t2.x()));
          if(minx - pt.x() > mindist)
            continue;
  
          // Do intersection
          if(rayXTriangleIntersect(pt, t0, t1, t2, n, intp)) {
            float dist = intp.x() - pt.x();
            if(dist > 0 and dist < mindist) {
              mindist = dist;
              minpt = intp;
              inside = (n.x() > 0 ? true : false);
            }
          }
        }
  
        // If no triangle found exit
        if(mindist == FLT_MAX)
          return;
        if(x == 0)
          pt = imageToWorld(imgBase + base + Point3i(x, y, z), step, shift);
        size_t xIdx = getOffset(x, y, z, size);
        while(pt.x() <= minpt.x() and x < size.x()) {
          if(inside)
            Dst[xIdx++] = 1;
          pt = imageToWorld(imgBase + base + Point3i(++x, y, z), step, shift);
        }
      } while(x < size.x());
    }
  };
  
  // Calculate mask contained by a volume
  struct nearMeshOP 
  {
    const Point3f step;
    const Point3f shift;
    const Point3f* pts;
    const Point3f* nrmls;
    const Point3i imgBase;
    const uint* gridIdx;
    const uint* gridCnt;
    const uint* gridPts;
    const Point3i gridSize;
    const Point3i gridPix;
    const float minDist;
    const float maxDist;
  
    Point3i size;
    Point3i base;
    Point3i stride;
    ushort* Dst;
  
    nearMeshOP(const Point3f& _step, const Point3f& _shift, const Point3f* _pts, const Point3f* _nrmls,
               const Point3i& _imgBase, const uint* _gridIdx, const uint* _gridCnt, const uint* _gridPts,
               const Point3i _gridSize, const Point3i& _gridPix, float _minDist, float _maxDist, DeviceVus& _Dst)
      : step(_step), shift(_shift), pts(_pts), nrmls(_nrmls), imgBase(_imgBase), gridIdx(_gridIdx),
        gridCnt(_gridCnt), gridPts(_gridPts), gridSize(_gridSize), gridPix(_gridPix), minDist(_minDist), 
        maxDist(_maxDist), Dst(devP(_Dst)) {}
  
    void setSizes(const Point3i& _base, const Point3i& _size)
    {
      base = _base;
      size = _size;
      stride = getStride(size);
    }
  
    __device__ __host__ void operator()(uint idx)
    {
      Point3f intp;
      int x, y, z;
      getXYZ(idx, size, x, y, z);
  
      // Calculate distance
      Point3f pt = imageToWorld(imgBase + base + Point3i(x, y, z), step, shift);
  
      // Grid information
      Point3i gridPos = divide(base + Point3i(x, y, z), gridPix);
      size_t idxG = getOffset(gridPos, gridSize);
  
      // Loop through triangle list to find closest triangle
      float dist = FLT_MAX;
      bool inside;
      int idxT = gridIdx[idxG];
      for(uint i = 0; i < gridCnt[idxG]; ++i, ++idxT) {
        Point3f p = pts[gridPts[idxT]];
        Point3f n = nrmls[gridPts[idxT]];
  
        float d = (p - pt).norm();
        if(d < dist) {
          dist = d;
          inside = ((p - pt) * n >= 0);
        }
      }
      // Simple case, single distance from mesh
      if(!inside)
        dist = -dist;
      if(dist >= minDist and dist <= maxDist)
        Dst[getOffset(x, y, z, size)] = 1;
    }
  };
  
  // Find min/max pixels for each slice
  struct pixelFindMinMaxXOP 
  {
    const Point3i base;
    const Point3i size;
    const Point3u imgSize;
    const Point3f step;
    const Point3f shift;
    const Point3f p;
    const Point3f pv;
    const Point2f width;
  
    Point2i* Dmin;
    Point2i* Dmax;
  
    pixelFindMinMaxXOP(const Point3i& _base, const Point3i& _size, const Point3u& _imgSize, const Point3f& _step,
                       const Point3f& _shift, const Point3f& _p, const Point3f& _pv, const Point2f& _width,
                       Point2i* _Dmin, Point2i* _Dmax)
      : base(_base), size(_size), imgSize(_imgSize), step(_step), shift(_shift), p(_p),
          pv(_pv), width(_width), Dmin(_Dmin), Dmax(_Dmax) {}
  
    __device__ __host__ void operator()(const uint x) const
    {
      Point3f pint;     // intersect point
      Point3f xp = imageToWorld(Point3i(base.x() + x, 0, 0), step, shift);
      float s;
  
      if(planeLineIntersect(xp, Point3f(1.0f, 0.0f, 0.0f), p, p + pv, s, pint)) {
        xp = pint;
        xp.z() -= width.x();
        xp.y() -= width.y();
        Point3i xi = worldToImage(xp, step, shift);
        Dmin[x].x() = trim(xi.z() - 1, base.z(), base.z() + size.z());
        Dmin[x].y() = trim(xi.y() - 1, base.y(), base.y() + size.y());
  
        xp = pint;
        xp.z() += width.x();
        xp.y() += width.y();
        xi = worldToImage(xp, step, shift);
        Dmax[x].x() = trim(xi.z() + 1, base.z(), base.z() + size.z());
        Dmax[x].y() = trim(xi.y() + 1, base.y(), base.y() + size.y());
      } else {
        Dmin[x] = Point2i(0, 0);
        Dmax[x] = Point2i(0, 0);
      }
    }
  };
  
  // Find min/max pixels for each slice
  struct pixelFindMinMaxYOP {
    const Point3i base;
    const Point3i size;
    const Point3u imgSize;
    const Point3f step;
    const Point3f shift;
    const Point3f p;
    const Point3f pv;
    const Point2f width;
  
    Point2i* Dmin;
    Point2i* Dmax;
  
    pixelFindMinMaxYOP(const Point3i& _base, const Point3i& _size, const Point3u& _imgSize, const Point3f& _step,
                       const Point3f& _shift, const Point3f& _p, const Point3f& _pv, const Point2f& _width,
                       Point2i* _Dmin, Point2i* _Dmax)
      : base(_base), size(_size), imgSize(_imgSize), step(_step), shift(_shift), p(_p), pv(_pv),
          width(_width), Dmin(_Dmin), Dmax(_Dmax) {}
  
    __device__ __host__ void operator()(const uint y) const
    {
      Point3f pint;     // intersect point
      Point3f yp = imageToWorld(Point3i(0, base.y() + y, 0), step, shift);
      float s;
  
      if(planeLineIntersect(yp, Point3f(0.0f, 1.0f, 0.0f), p, p + pv, s, pint)) {
        yp = pint;
        yp.x() -= width.x();
        yp.z() -= width.y();
        Point3i yi = worldToImage(yp, step, shift);
        Dmin[y].x() = trim(yi.x() - 1, base.x(), base.x() + size.x());
        Dmin[y].y() = trim(yi.z() - 1, base.z(), base.z() + size.z());
  
        yp = pint;
        yp.x() += width.x();
        yp.z() += width.y();
        yi = worldToImage(yp, step, shift);
        Dmax[y].x() = trim(yi.x() + 1, base.x(), base.x() + size.x());
        Dmax[y].y() = trim(yi.z() + 1, base.z(), base.z() + size.z());
      } else {
        Dmin[y] = Point2i(0, 0);
        Dmax[y] = Point2i(0, 0);
      }
    }
  };
  
  // Find min/max pixels for each slice
  struct pixelFindMinMaxZOP 
  {
    const Point3i base;
    const Point3i size;
    const Point3u imgSize;
    const Point3f step;
    const Point3f shift;
    const Point3f p;
    const Point3f pv;
    const Point2f width;
  
    Point2i* Dmin;
    Point2i* Dmax;
  
    pixelFindMinMaxZOP(const Point3i& _base, const Point3i& _size, const Point3u& _imgSize, const Point3f& _step,
                       const Point3f& _shift, const Point3f& _p, const Point3f& _pv, const Point2f& _width,
                       Point2i* _Dmin, Point2i* _Dmax)
      : base(_base), size(_size), imgSize(_imgSize), step(_step), shift(_shift), p(_p),
          pv(_pv), width(_width), Dmin(_Dmin), Dmax(_Dmax) {}

    __device__ __host__ void operator()(const uint z) const
    {
      Point3f pint;     // intersect point
      Point3f zp = imageToWorld(Point3i(0, 0, base.z() + z), step, shift);
      float s;
  
      if(planeLineIntersect(zp, Point3f(0.0f, 0.0f, 1.0f), p, p + pv, s, pint)) {
        zp = pint;
        zp.x() -= width.x();
        zp.y() -= width.y();
        Point3i zi = worldToImage(zp, step, shift);
        Dmin[z].x() = trim(zi.x() - 1, base.x(), base.x() + size.x());
        Dmin[z].y() = trim(zi.y() - 1, base.y(), base.y() + size.y());
  
        zp = pint;
        zp.x() += width.x();
        zp.y() += width.y();
        zi = worldToImage(zp, step, shift);
        Dmax[z].x() = trim(zi.x() + 1, base.x(), base.x() + size.x());
        Dmax[z].y() = trim(zi.y() + 1, base.y(), base.y() + size.y());
  
      } else {
        Dmin[z] = Point2i(0, 0);
        Dmax[z] = Point2i(imgSize.x(), imgSize.y());
      }
    }
  };
  
  // Find min/max pixels for each slice
  struct voxelEditXOP {
    const uint dbase;
    const uint gbase;
    const Point3u size;
    const Point3u imgSize;
    const Point3f step;
    const Point3f shift;
    const Point3f p;
    const Point3f pv;
    const float radius;
    const ushort color;
    const Point3u clipDo;
    const Point4f* pn;
  
    const Point2i* Dmin;
    const Point2i* Dmax;
    ushort* Dpix;
  
    voxelEditXOP(uint _dbase, uint _gbase, const Point3u& _size, const Point3u& _imgSize, const Point3f& _step,
                 const Point3f& _shift, const Point3f& _p, const Point3f& _pv, float _radius, ushort _color,
                 const Point3u& _clipDo, const Point4f* _pn, const Point2i* _Dmin, const Point2i* _Dmax, ushort* _Dpix)
      : dbase(_dbase), gbase(_gbase), size(_size), imgSize(_imgSize), step(_step), shift(_shift), p(_p), pv(_pv), 
          radius(_radius), color(_color), clipDo(_clipDo), pn(_pn), Dmin(_Dmin), Dmax(_Dmax), Dpix(_Dpix) {}
  
    __device__ __host__ void operator()(const uint idx) const
    {
      uint yb = idx % size.y();
      uint zb = idx / size.y();
      for(int x = 0; x < (int)size.x(); ++x) {
        uint dx = x + dbase;
        uint z = zb + Dmin[dx].x();
        uint y = yb + Dmin[dx].y();
        if((int)z >= Dmax[dx].x() || (int)y >= Dmax[dx].y())
          continue;
  
        bool erase = true;
        Point3f wp = imageToWorld(Point3i(gbase + dx, y, z), step, shift);
        if(clipPoint(wp, clipDo, pn))
          erase = false;
  
        // Test if within radius
        if(distLinePoint(p - pv, p + pv, wp, false) >= radius)
          erase = false;
  
        // Save result
        if(erase)
          Dpix[offset(yb, zb, x, size.y(), size.z())] = color;
      }
    }
  };
  
  // Find min/max pixels for each slice
  struct voxelEditYOP 
  {
    const uint dbase;
    const uint gbase;
    const Point3u size;
    const Point3u imgSize;
    const Point3f step;
    const Point3f shift;
    const Point3f p;
    const Point3f pv;
    const float radius;
    const ushort color;
    const Point3u clipDo;
    const Point4f* pn;
  
    const Point2i* Dmin;
    const Point2i* Dmax;
    ushort* Dpix;
  
    voxelEditYOP(uint _dbase, uint _gbase, const Point3u& _size, const Point3u& _imgSize, const Point3f& _step,
                 const Point3f& _shift, const Point3f& _p, const Point3f& _pv, float _radius, ushort _color,
                 const Point3u& _clipDo, const Point4f* _pn, const Point2i* _Dmin, const Point2i* _Dmax, ushort* _Dpix)
      : dbase(_dbase), gbase(_gbase), size(_size), imgSize(_imgSize), step(_step), shift(_shift), p(_p), pv(_pv),
          radius(_radius), color(_color), clipDo(_clipDo), pn(_pn), Dmin(_Dmin), Dmax(_Dmax), Dpix(_Dpix) {}
  
    __device__ __host__ void operator()(const uint idx) const
    {
      uint xb = idx % size.x();
      uint zb = idx / size.x();
      for(int y = 0; y < (int)size.y(); y++) {
        uint dy = dbase + y;
        uint x = xb + Dmin[dy].x();
        uint z = zb + Dmin[dy].y();
        if((int)x >= Dmax[dy].x() || (int)z >= Dmax[dy].y())
          continue;
  
        bool erase = true;
        Point3f wp = imageToWorld(Point3i(x, gbase + dy, z), step, shift);
        if(clipPoint(wp, clipDo, pn))
          erase = false;
  
        // Test if within radius
        if(distLinePoint(p - pv, p + pv, wp, false) >= radius)
          erase = false;
  
        // Save result
        if(erase)
          Dpix[offset(xb, zb, y, size.x(), size.z())] = color;
      }
    }
  };
  
  // Find min/max pixels for each slice
  struct voxelEditZOP 
  {
    const uint dbase;
    const uint gbase;
    const Point3u size;
    const Point3u imgSize;
    const Point3f step;
    const Point3f shift;
    const Point3f p;
    const Point3f pv;
    const float radius;
    const ushort color;
    const Point3u clipDo;
    const Point4f* pn;
  
    Point2i* Dmin;
    Point2i* Dmax;
    ushort* Dpix;
  
    voxelEditZOP(uint _dbase, uint _gbase, const Point3u& _size, const Point3u& _imgSize, const Point3f& _step,
                 const Point3f& _shift, const Point3f& _p, const Point3f& _pv, float _radius, ushort _color,
                 const Point3u& _clipDo, const Point4f* _pn, Point2i* _Dmin, Point2i* _Dmax, ushort* _Dpix)
      : dbase(_dbase), gbase(_gbase), size(_size), imgSize(_imgSize), step(_step), 
        shift(_shift), p(_p), pv(_pv), radius(_radius), color(_color), clipDo(_clipDo), 
        pn(_pn), Dmin(_Dmin), Dmax(_Dmax), Dpix(_Dpix) {}
  
    __device__ __host__ void operator()(uint idx) const
    {
      uint xb = idx % size.x();
      uint yb = idx / size.x();
      for(int z = 0; z < (int)size.z(); z++) {
        uint dz = z + dbase;
        uint x = xb + Dmin[dz].x();
        uint y = yb + Dmin[dz].y();
        if((int)x >= Dmax[dz].x() || (int)y >= Dmax[dz].y())
          continue;
  
        bool erase = true;
        // Test against clip planes
        Point3f wp = imageToWorld(Point3i(x, y, gbase + dz), step, shift);
        if(clipPoint(wp, clipDo, pn))
          erase = false;
  
        // Test if within radius
        if(distLinePoint(p - pv, p + pv, wp, false) >= radius)
          erase = false;
  
        // Save result
        if(erase)
          Dpix[offset(xb, yb, z, size.x(), size.y())] = color;
      }
    }
  };

  // ---------- Externally exported functions ----------
  
  // Average voxels (separable)
  cuda_EXPORT
  int averageGPU(const Point3i& imgSize, const Point3i& radius, const HostVus& srcdata, HostVus& dstdata)
  {
    HoldMem hm;
  
    AverageKernel op = AverageKernel();
    return processSepOP(op, op, op, imgSize, radius, srcdata, dstdata);
  }
  
  // Median of voxels (non-separable)
  cuda_EXPORT
  int medianGPU(const Point3i& imgSize, const Point3i& radius, const HostVus& srcdata, HostVus& dstdata)
  {
    HoldMem hm;
  
    MedianKernel op = MedianKernel();
    return processNonSepOP(op, imgSize, radius, srcdata, dstdata);
  }
  
  // Median of medians (separable)
  cuda_EXPORT
  int medianOfMediansGPU(const Point3i& imgSize, const Point3i& radius, const HostVus& srcdata, HostVus& dstdata)
  {
    HoldMem hm;
  
    MedianKernel op = MedianKernel();
    return processSepOP(op, op, op, imgSize, radius, srcdata, dstdata);
  }
  
  // Dilate voxels (separable)
  cuda_EXPORT
  int dilateGPU(const Point3i& imgSize, const Point3i& radius, bool roundNhbd, const HostVus& srcdata, HostVus& dstdata)
  {
    HoldMem hm;
  
    DeviceVf dMask;
  
    if(roundNhbd) {
      HostVf mask;
      getRoundMask3D(radius, mask);
      DeviceVf dMask(mask);
      return processNonSepOP(DilateMaskKernel(devP(dMask)), imgSize, radius, srcdata, dstdata);
    } else {
      DilateKernel op = DilateKernel();
      return processSepOP(op, op, op, imgSize, radius, srcdata, dstdata);
    }
  }
  
  // Erode voxels (separable)
  cuda_EXPORT
  int erodeGPU(const Point3i& imgSize, const Point3i& radius, bool byLabel, bool roundNhbd, 
  								                                          const HostVus& srcdata, HostVus& dstdata)
  {
    HoldMem hm;
  
    DeviceVf dMask;
  
    if(byLabel) {
      ErodeLabelsKernel op = ErodeLabelsKernel();
      return processSepOP(op, op, op, imgSize, radius, srcdata, dstdata);
    } else if(roundNhbd) { // So far not by label
      HostVf mask;
      getRoundMask3D(radius, mask);
      DeviceVf dMask(mask);
      return processNonSepOP(ErodeMaskKernel(devP(dMask)), imgSize, radius, srcdata, dstdata);
    } else {
      ErodeKernel op = ErodeKernel();
      return processSepOP(op, op, op, imgSize, radius, srcdata, dstdata);
    }
  }
  
  // Guassian blur voxels (separable)
  cuda_EXPORT
  int gaussianBlurGPU(const Point3i& imgSize, const Point3f& imgStep, const Point3f& sigma, const HostVus& srcdata,
                      HostVus& dstdata)
  {
    HoldMem hm;
  
    Point3i radius(0, 0, 0);
    DeviceVf dMask[3];
  
    // Setup masks
    for(int i = 0; i < 3; ++i) {
      HostVf mask;
      getGaussianMask1D(imgStep[i], sigma[i], mask, radius[i]);
      normalizeMaskScale(mask);
      dMask[i] = mask;
    }
    return processSepOP(MaskKernel<ushort>(devP(dMask[0])), MaskKernel<ushort>(devP(dMask[1])), MaskKernel<ushort>(devP(dMask[2])), imgSize,
                        radius, srcdata, dstdata);
  }
  
  // Sharpen voxels (separable)
  cuda_EXPORT
  int sharpenGPU(const Point3i& imgSize, const Point3f& imgStep, const Point3f& sigma, const float amount,
                 const HostVus& srcdata, HostVus& dstdata)
  {
    HoldMem hm;
  
    Point3i radius(0, 0, 0);
    DeviceVf dMask[3];
  
    // Setup masks
    for(int i = 0; i < 3; ++i) {
      HostVf mask;
      getGaussianMask1D(imgStep[i], sigma[i], mask, radius[i]);
      normalizeMaskScale(mask);
  
      if(radius[i] > 0) {
        // weight center, and flip sign
        float normalize = 0;
        for(int j = 0; j < (int)mask.size(); ++j)
          if(j != radius[i]) {
            normalize += mask[j];
            mask[j] *= -amount;
            normalize -= mask[j];
          }
        mask[radius[i]] += normalize;
      }
      dMask[i] = mask;
    }
    return processSepOP(MaskKernel<ushort>(devP(dMask[0])), MaskKernel<ushort>(devP(dMask[1])), MaskKernel<ushort>(devP(dMask[2])), imgSize,
                        radius, srcdata, dstdata);
  }
  
  // Apply a separable kernel
  cuda_EXPORT
  int applyKernelGPU(const Point3i& imgSize, const HostVf& kernelX, const HostVf& kernelY, const HostVf& kernelZ,
                     const HostVus& srcdata, HostVus& dstdata)
  {
    HoldMem hm;
  
    Point3i radius((kernelX.size() - 1) / 2, (kernelY.size() - 1) / 2, (kernelZ.size() - 1) / 2);
    DeviceVf dKernelX(kernelX);
    DeviceVf dKernelY(kernelY);
    DeviceVf dKernelZ(kernelZ);
  
    return processSepOP(MaskKernel<ushort>(devP(dKernelX)), MaskKernel<ushort>(devP(dKernelY)), MaskKernel<ushort>(devP(dKernelZ)), imgSize,
                        radius, srcdata, dstdata);
  }
  
  // Color gradient (in Z only)
  cuda_EXPORT
  int colorGradGPU(const Point3i& imgSize, float div, const HostVus& srcData, HostVus& dstData)
  {
    HoldMem hm;
  
    HostVus Hsrc(srcData);
    DeviceVus Dsrc, Ddst;
    processOP proc(imgSize, 1, sizeof(ushort) * 2);
    do {
      proc.write(Hsrc, Dsrc);
      proc.alloc(Ddst);
      if(proc.launch(colorGradOP(div, &Dsrc, &Ddst), XY))
        return 1;
      proc.read(Ddst, dstData);
    } while(proc.next());
  
    return 0;
  }
  
  // Normalized edge detect
  cuda_EXPORT
  int edgeDetectGPU(const Point3i& imgSize, ushort lowthresh, ushort highthresh, float mult, 
    ushort fillval, const HostVus& srcData, HostVus& dstData)
  {
    HoldMem hm;
  
    HostVus Hsrc(srcData.size()), Hdst(srcData.size());
    // Flip Y and Z axes for processing
    Point3i flipImgSize;
    flipImage(YZ, imgSize, flipImgSize, srcData, Hsrc);
  
    DeviceVus Dsrc, Ddst;
    processOP proc(flipImgSize, 0, sizeof(ushort) * 2);
    do {
      proc.write(Hsrc, Dsrc);
      proc.alloc(Ddst);
      if(proc.launch(edgeDetectOP(lowthresh, highthresh, mult, fillval, &Dsrc, &Ddst), XZ))
        return 1;
      proc.read(Ddst, Hdst);
    } while(proc.next());
  
    flipImage(YZ, flipImgSize, flipImgSize, Hdst, dstData);
  
    return 0;
  }
  
  // Clip stack to clipping planes
  cuda_EXPORT
  int clipStackGPU(const Point3i& imgSize, const Point3f& step, const Point3f& shift, const Point3u& clipDo, HostV4f& pn,
                   const HostVus& srcData, HostVus& dstData)
  {
    HoldMem hm;
  
    // If all clip planes disabled, just copy data
    if(clipDo.x() == 0 && clipDo.y() == 0 && clipDo.z() == 0) {
      dstData = srcData;
      return (0);
    }
    dstData.resize(srcData.size());
    DeviceV4f Dpn(pn);   // Put clip planes on device
    DeviceVus Dsrc, Ddst;
    processOP proc(imgSize, 0, sizeof(ushort) * 2);
    do {
      proc.write(srcData, Dsrc);
      proc.alloc(Ddst);
      if(proc.launch(clipStackOP(step, shift, clipDo, Dpn, &Dsrc, &Ddst), XY))
        return 1;
      proc.read(Ddst, dstData);
    } while(proc.next());
  
    return 0;
  }
  
  // Find out if pixels are inside mesh (closed), return a list of pixels inside
  cuda_EXPORT
  int insideMeshGPU(const Point3i& base, const Point3i& size, const Point3f& step, 
    const Point3f& shift, const HostV3f& pts, const HostV3u& tris, const HostV3f& nrmls, HostVus& dstData)
  {
    HoldMem hm;
  
    if(size_t(size.x()) * size.y() * size.z() <= 0)
      return errMsg("insideMeshGPU:Empty region");
  
    // Create YZ grid
    const Point3i gridPix(1, 32, 8);
    const Point3i gridSize = divide(size, gridPix) + Point3i(1, 1, 1);
    std::vector<std::set<int> > grid;
    grid.resize(gridSize.y() * gridSize.z());
  
    for(size_t i = 0; i < tris.size(); i++) {
      Point3f t0 = pts[tris[i].x()];
      Point3f t1 = pts[tris[i].y()];
      Point3f t2 = pts[tris[i].z()];
  
      // Triangle min max in world coordinates
      Point3f minW = min(t0, min(t1, t2));
      Point3f maxW = max(t0, max(t1, t2));
      // Min max in image coords
      Point3i minI = worldToImage(minW, step, shift) - base;
      Point3i maxI = worldToImage(maxW, step, shift) - base;
      Point3i minG = divide(minI, gridPix);
      Point3i maxG = divide(maxI, gridPix);
      // Triangle is outside grid
      if(maxG.y() < 0 or minG.y() >= gridSize.y())
        continue;
      if(maxG.z() < 0 or minG.z() >= gridSize.z())
        continue;
      minG = trim(minG, Point3i(0, 0, 0), gridSize - Point3i(1, 1, 1));
      maxG = trim(maxG, Point3i(0, 0, 0), gridSize - Point3i(1, 1, 1));
      // Add triangle to all grid squares
      for(int z = minG.z(); z <= maxG.z(); z++)
        for(int y = minG.y(); y <= maxG.y(); y++)
          grid[z * gridSize.y() + y].insert(i);
    }
    int maxtris = 0;
    int tottris = 0;
    for(size_t i = 0; i < grid.size(); i++) {
      int numtris = grid[i].size();
      tottris += numtris;
      if(maxtris < numtris)
        maxtris = numtris;
    }
    cout << "Grid Size:" << gridSize.y() << "x" << gridSize.z() << " MaxTris:" << maxtris << " TotTris:" << tottris
         << " Tris:" << tris.size() << endl;
    HostVu hGridIdx(grid.size());
    HostVu hGridCnt(grid.size());
    HostVu hGridTris(tottris);
    int idx = 0;
    for(size_t i = 0; i < grid.size(); i++) {
      hGridIdx[i] = idx;
      hGridCnt[i] = grid[i].size();
      std::copy(grid[i].begin(), grid[i].end(), hGridTris.begin() + idx);
      idx += grid[i].size();
    }
  
    // Check memory, we need room for points, triangles, normals, grid, and at least one slice
    size_t reqmem = (nrmls.size() + pts.size()) * sizeof(Point3f) + tris.size() * sizeof(Point3u)
      + (hGridIdx.size() + hGridCnt.size() + hGridTris.size()) * sizeof(uint);
  
    if(reqmem + 2 * size.x() * size.y() >= userMem())
      throw(string("InsideMeshGPU::Not enough GPU memory to store triangle list, try increasing CudaHoldMem"));
  
    // Put points, triangles, normals, and grid on device
    DeviceV3f dPts(pts);
    DeviceV3u dTris(tris);
    DeviceV3f dNrmls(nrmls);
    DeviceVu dGridIdx(hGridIdx);
    DeviceVu dGridCnt(hGridCnt);
    DeviceVu dGridTris(hGridTris);
    DeviceVus Ddst;
    checkCudaError("InsideMeshGPU::Allocating device buffers for triangle list.");
  
    processOP proc(size, 0, sizeof(ushort));
    do {
      proc.alloc(Ddst);
      if(proc.launch(insideMeshOP(step, shift, devP(dPts), devP(dTris), devP(dNrmls), 
        base, devP(dGridIdx), devP(dGridCnt), devP(dGridTris), gridSize, gridPix, Ddst),
                     YZ))
        return 1;
      proc.read(Ddst, dstData);
    } while(proc.next());
  
    return 0;
  }
  
  // Find out if pixels are inside mesh (closed), and return a mask
  cuda_EXPORT
  int nearMeshGPU(const Point3i& base, const Point3i& size, const Point3f& step, 
    const Point3f& shift, const HostV3f& pts, const HostV3f& nrmls, float minDist, 
    float maxDist, HostVus& dstData)
  {
    HoldMem hm;
  
    if(size_t(size.x()) * size.y() * size.z() <= 0)
      return errMsg("nearMeshGPU:Empty region");
    if(minDist >= maxDist)
      return errMsg("nearMeshGPU:Min Distance must be less than Max Distance");
  
    // Create XYZ grid
    float distance = max(fabs(minDist), fabs(maxDist));
    Point3f distance3f(distance, distance, distance);
    const Point3i gridPix = max(Point3i(32, 32, 8), Point3i(divide(distance3f, step)));
    const Point3i gridSize = divide(size, gridPix) + Point3i(1, 1, 1);
    std::vector<std::set<int> > grid;
    grid.resize(gridSize.x() * gridSize.y() * gridSize.z());
  
    for(size_t i = 0; i < pts.size(); i++) {
      Point3f p = pts[i];
  
      // Point min max in world coordinates
      Point3f minW = p - distance3f;
      Point3f maxW = p + distance3f;
      // Min max in image coords
      Point3i minI = worldToImage(minW, step, shift) - base;
      Point3i maxI = worldToImage(maxW, step, shift) - base;
      // Find grid coordinates
      Point3i minG = divide(minI, gridPix);
      Point3i maxG = divide(maxI, gridPix);
      // Point influence is outside grid
      if(maxG.x() < 0 or minG.x() >= gridSize.x())
        continue;
      if(maxG.y() < 0 or minG.y() >= gridSize.y())
        continue;
      if(maxG.z() < 0 or minG.z() >= gridSize.z())
        continue;
      minG = trim(minG, Point3i(0, 0, 0), gridSize - Point3i(1, 1, 1));
      maxG = trim(maxG, Point3i(0, 0, 0), gridSize - Point3i(1, 1, 1));
      // Add point to all grid squares
      for(int z = minG.z(); z <= maxG.z(); z++)
        for(int y = minG.y(); y <= maxG.y(); y++)
          for(int x = minG.x(); x <= maxG.x(); x++)
            grid[getOffset(x, y, z, gridSize)].insert(i);
    }
    int maxpts = 0;
    int totpts = 0;
    for(size_t i = 0; i < grid.size(); i++) {
      int numpts = grid[i].size();
      totpts += numpts;
      if(maxpts < numpts)
        maxpts = numpts;
    }
    cout << "Grid Size:" << gridSize.x() << "x" << gridSize.y() << "x" << gridSize.z() << " MaxPts:" << maxpts
         << " TotPts:" << totpts << " Pts:" << pts.size() << endl;
    HostVu hGridIdx(grid.size());
    HostVu hGridCnt(grid.size());
    HostVu hGridPts(totpts);
    int idx = 0;
    for(size_t i = 0; i < grid.size(); i++) {
      hGridIdx[i] = idx;
      hGridCnt[i] = grid[i].size();
      std::copy(grid[i].begin(), grid[i].end(), hGridPts.begin() + idx);
      idx += grid[i].size();
    }
  
    // Check memory, we need room for points, triangles, normals, grid, and at least one slice
    size_t reqmem = (nrmls.size() + pts.size()) * sizeof(Point3f)
      + (hGridIdx.size() + hGridCnt.size() + hGridPts.size()) * sizeof(uint);
  
    if(reqmem + 2 * size.x() * size.y() >= userMem())
      throw(string("NearMeshGPU::Not enough GPU memory to store point list, try increasing CudaHoldMem"));
  
    // Put points, triangles, normals, and grid on device
    DeviceV3f dPts(pts);
    DeviceV3f dNrmls(nrmls);
    DeviceVu dGridIdx(hGridIdx);
    DeviceVu dGridCnt(hGridCnt);
    DeviceVu dGridPts(hGridPts);
    DeviceVus Ddst;
    checkCudaError("NearMeshGPU::Allocating device buffers for point list");
  
    processOP proc(size, 0, sizeof(ushort));
    do {
      proc.alloc(Ddst);
      if(proc.launch(nearMeshOP(step, shift, devP(dPts), devP(dNrmls), base, devP(dGridIdx), devP(dGridCnt),
                                devP(dGridPts), gridSize, gridPix, minDist, maxDist, Ddst), XYZ))
        return 1;
      proc.read(Ddst, dstData);
    } while(proc.next());
  
    return 0;
  }
  
  // Edit pixel data in 3D, erase, fill, and seed
  cuda_EXPORT
  int voxelEditGPU(const Point3i& base, const Point3i& size, const Point3i& imgSize, const Point3f& step,
                   const Point3f& shift, const Point3f& p, const Point3f& pv, float radius, ushort color,
                   const Point3u& clipDo, const HostV4f& pn, HostVus& data)
  {
    HoldMem hm;
  
    bool err = 0;
  
    // Check bounds
    if(checkBounds(imgSize, base, size))
      throw(string("Bounds error (or empty)"));
    else
      try {
        // Max image coordinate
        size_t imgMax = size_t(imgSize.x()) * imgSize.y() * imgSize.z();
  
        // Declare device vectors
        DeviceV4f Dpn(pn);       // put clip planes on device
        checkCudaError("Allocating buffer for clipping planes");
  
        // Favor Z direction
        float vxy = Point2f(pv.x(), pv.y()).norm();
        if(vxy < .1 || fabs(pv.z() / vxy) > .4) {       // Process in z direction
          // Calculate memory
          size_t mem = memLeft(size.z() * 2 * sizeof(Point2i), 4 * MByte, userMem());
  
          // Allocate slice min/max vectors
          DeviceV2i Dmin(size.z()), Dmax(size.z());
          checkCudaError("Allocating min/max device vectors");
  
          // Cross product with plane projected view vector (ie, onto xz && yz)
          Point3f px(pv.x(), 0, pv.z());
          px.normalize();
          Point3f py(0, pv.y(), pv.z());
          py.normalize();
          Point2f width(radius / sqrt(px.z() * px.z()), radius / sqrt(py.z() * py.z()));
  
          // Launch threads on section
          thrust::for_each(DCountIter(0), DCountIter(size.z()),
                           pixelFindMinMaxZOP(base, size, Point3u(imgSize), step, shift, p, pv, width, devP(Dmin),
                                              devP(Dmax)));
          checkCudaError("pixelFindMinMaxZOP launch");
  
          // Get regions to process
          HostV2i Hmin(Dmin), Hmax(Dmax);
          checkCudaError("Min/max copy to host");
  
          // Find max size && which slices to process
          Point3i sz(0, 0, 0);
          for(int zIdx = 0; zIdx < size.z(); zIdx++) {
            int xsz = Hmax[zIdx].x() - Hmin[zIdx].x();
            int ysz = Hmax[zIdx].y() - Hmin[zIdx].y();
  
            if(xsz > sz.x())
              sz.x() = xsz;
            if(ysz > sz.y())
              sz.y() = ysz;
          }
          if(sz.x() > 0 && sz.y() > 0) {
            // Collect pixel data
            HostVus Hpix(size_t(sz.x()) * sz.y() * size.z());
            size_t idx = 0;
            for(int z = 0; z < size.z(); z++)
              for(int y = Hmin[z].y(); y < Hmin[z].y() + sz.y(); y++)
                for(int x = Hmin[z].x(); x < Hmin[z].x() + sz.x(); x++) {
                  if(x < Hmax[z].x() && y < Hmax[z].y()) {
                    size_t gidx = offset(x, y, z + base.z(), imgSize.x(), imgSize.y());
                    if(gidx < imgMax)
                      Hpix[idx] = data[gidx];
                  }
                  idx++;
                }
  
            // Process pixels
            size_t zIdx = 0;
            while((sz.z() = getStep(zIdx, size.z(), sz.x() * sz.y(), sizeof(ushort), mem)) > 0) {
              // Put data on device
              DeviceVus Dpix(Hpix.begin() + zIdx * sz.x() * sz.y(),
                          Hpix.begin() + size_t(zIdx + sz.z()) * sz.x() * sz.y());
              checkCudaError("Copy pixels to device");
  
              // Process pixels, || far only erase
              thrust::for_each(DCountIter(0), DCountIter(sz.x() * sz.y()),
                               voxelEditZOP(zIdx, base.z(), Point3u(sz), Point3u(imgSize), step, shift, p, pv,
                                            radius, color, clipDo, devP(Dpn), devP(Dmin), devP(Dmax),
                                            devP(Dpix)));
              checkCudaError("pixEraseZOP launch");
  
              // Copy data back to host
              thrust::copy(Dpix.begin(), Dpix.end(), Hpix.begin() + size_t(zIdx) * sz.x() * sz.y());
              checkCudaError("Copying pixels to host");
  
              zIdx += sz.z();
            }
  
            // Update results
            idx = 0;
            for(int z = 0; z < size.z(); z++)
              for(int y = Hmin[z].y(); y < Hmin[z].y() + sz.y(); y++)
                for(int x = Hmin[z].x(); x < Hmin[z].x() + sz.x(); x++) {
                  if(x < Hmax[z].x() && y < Hmax[z].y()) {
                    size_t gidx = offset(x, y, base.z() + z, imgSize.x(), imgSize.y());
                    if(gidx < imgMax)
                      data[gidx] = Hpix[idx];
                  }
                  idx++;
                }
          }
        } else if(fabs(pv.y()) > fabs(pv.x())) {       // Process Y direction
          // Calculate memory
          size_t mem = memLeft(size.y() * 2 * sizeof(Point2i), 4 * MByte, userMem());
  
          // Allocate slice min/max vectors
          DeviceV2i Dmin(size.y()), Dmax(size.y());
          checkCudaError("Allocating min/max device vectors");
  
          // Cross product with plane projected view vector (ie, onto xy && zy)
          Point3f px = Point3f(pv.x(), pv.y(), 0).normalize();
          Point3f pz = Point3f(0, pv.y(), pv.z()).normalize();
          Point2f width = Point2f(radius / sqrt(px.y() * px.y()), radius / sqrt(pz.y() * pz.y()));
  
          // Launch threads on section
          thrust::for_each(DCountIter(0), DCountIter(size.y()),
                           pixelFindMinMaxYOP(base, size, Point3u(imgSize), step, shift, p, pv, width, devP(Dmin),
                                              devP(Dmax)));
          checkCudaError("pixelFindMinMaxYOP launch");
  
          // Get regions to process
          HostV2i Hmin(Dmin), Hmax(Dmax);
          checkCudaError("Min/max copy to host");
  
          // Find max size && which slices to process
          Point3i sz(0, 0, 0);
          for(int yIdx = 0; yIdx < size.y(); yIdx++) {
            int xsz = Hmax[yIdx].x() - Hmin[yIdx].x();
            int zsz = Hmax[yIdx].y() - Hmin[yIdx].y();
  
            if(xsz > sz.x())
              sz.x() = xsz;
            if(zsz > sz.z())
              sz.z() = zsz;
          }
          if(sz.x() > 0 && sz.z() > 0) {
            // Collect pixel data
            HostVus Hpix(size_t(sz.x()) * size.y() * sz.z());
            size_t idx = 0;
            for(int y = 0; y < size.y(); y++)
              for(int z = Hmin[y].y(); z < Hmin[y].y() + sz.z(); z++)
                for(int x = Hmin[y].x(); x < Hmin[y].x() + sz.x(); x++) {
                  if(x < Hmax[y].x() && z < Hmax[y].y()) {
                    size_t gidx = offset(x, base.y() + y, z, imgSize.x(), imgSize.y());
                    if(gidx < imgMax)
                      Hpix[idx] = data[gidx];
                  }
                  idx++;
                }
  
            // Process pixels
            int yIdx = 0;
            while((sz.y() = getStep(yIdx, size.y(), sz.x() * sz.z(), sizeof(ushort), mem)) > 0) {
              // Put data on device
              DeviceVus Dpix(Hpix.begin() + size_t(yIdx) * sz.x() * sz.z(),
                          Hpix.begin() + size_t(yIdx + sz.y()) * sz.x() * sz.z());
              checkCudaError("Copy pixels to device");
  
              // Process pixels
              thrust::for_each(DCountIter(0), DCountIter(sz.x() * sz.z()),
                               voxelEditYOP(yIdx, base.y(), Point3u(sz), Point3u(imgSize), step, shift, p, pv,
                                            radius, color, clipDo, devP(Dpn), devP(Dmin), devP(Dmax),
                                            devP(Dpix)));
              checkCudaError("pixEraseYOP launch");
  
              // Copy data back to host
              thrust::copy(Dpix.begin(), Dpix.end(), Hpix.begin() + size_t(yIdx) * sz.x() * sz.z());
              checkCudaError("Copying pixels to host");
  
              yIdx += sz.y();
            }
  
            // Update reuslts
            idx = 0;
            for(int y = 0; y < size.y(); y++)
              for(int z = Hmin[y].y(); z < Hmin[y].y() + sz.z(); z++)
                for(int x = Hmin[y].x(); x < Hmin[y].x() + sz.x(); x++) {
                  if(x < Hmax[y].x() && z < Hmax[y].y()) {
                    size_t gidx = offset(x, base.y() + y, z, imgSize.x(), imgSize.y());
                    if(gidx < imgMax)
                      data[gidx] = Hpix[idx];
                  }
                  idx++;
                }
          }
        } else {       // else do X direction
          // Calculate memory
          size_t mem = memLeft(size.x() * 2 * sizeof(Point2i), 4 * MByte, userMem());
  
          // Allocate slice min/max vectors
          DeviceV2i Dmin(size.x()), Dmax(size.x());
          checkCudaError("Allocating min/max device vectors");
  
          // Cross product with plane projected view vector (i,.e, onto xz && yz)
          Point3f pz = Point3f(pv.x(), 0, pv.z()).normalize();
          Point3f py = Point3f(pv.x(), pv.y(), 0).normalize();
          Point2f width = Point2f(radius / sqrt(pz.x() * pz.x()), radius / sqrt(py.x() * py.x()));
  
          // Launch threads on section
          thrust::for_each(DCountIter(0), DCountIter(size.x()),
                           pixelFindMinMaxXOP(base, size, Point3u(imgSize), step, shift, p, pv, width, devP(Dmin),
                                              devP(Dmax)));
          checkCudaError("pixelFindMinMaxXOP launch");
  
          // Get regions to process
          HostV2i Hmin(Dmin), Hmax(Dmax);
          checkCudaError("Min/max copy to host");
  
          // Find max size && which slices to process
          Point3i sz(0, 0, 0);
          for(int xIdx = 0; xIdx < size.x(); xIdx++) {
            int zsz = Hmax[xIdx].x() - Hmin[xIdx].x();
            int ysz = Hmax[xIdx].y() - Hmin[xIdx].y();
  
            if(zsz > sz.z())
              sz.z() = zsz;
            if(ysz > sz.y())
              sz.y() = ysz;
          }
          if(sz.y() > 0 && sz.z() > 0) {
            // Collect pixel data
            HostVus Hpix(size_t(size.x()) * sz.y() * sz.z());
            size_t idx = 0;
            for(int x = 0; x < size.x(); x++)
              for(int z = Hmin[x].x(); z < Hmin[x].x() + sz.z(); z++)
                for(int y = Hmin[x].y(); y < Hmin[x].y() + sz.y(); y++) {
                  if(z < Hmax[x].x() && y < Hmax[x].y()) {
                    size_t gidx = offset(base.x() + x, y, z, imgSize.x(), imgSize.y());
                    if(gidx < imgMax)
                      Hpix[idx] = data[gidx];
                  }
                  idx++;
                }
  
            // Process pixels
            int xIdx = 0;
            while((sz.x() = getStep(xIdx, size.x(), sz.y() * sz.z(), sizeof(ushort), mem)) > 0) {
              // Put data on device
              DeviceVus Dpix(Hpix.begin() + xIdx * sz.y() * sz.z(),
                          Hpix.begin() + (xIdx + sz.x()) * sz.y() * sz.z());
              checkCudaError("Copy pixels to device");
  
              // Process pixels
              thrust::for_each(DCountIter(0), DCountIter(sz.y() * sz.z()),
                               voxelEditXOP(xIdx, base.x(), Point3u(sz), Point3u(imgSize), step, shift, p, pv,
                                            radius, color, clipDo, devP(Dpn), devP(Dmin), devP(Dmax),
                                            devP(Dpix)));
              checkCudaError("pixEraseXOP launch");
  
              // Copy data back to host
              thrust::copy(Dpix.begin(), Dpix.end(), Hpix.begin() + xIdx * sz.y() * sz.z());
              checkCudaError("Copying pixels to host");
  
              xIdx += sz.x();
            }
  
            // Update reuslts
            idx = 0;
            for(int x = 0; x < size.x(); x++)
              for(int z = Hmin[x].x(); z < Hmin[x].x() + sz.z(); z++)
                for(int y = Hmin[x].y(); y < Hmin[x].y() + sz.y(); y++) {
                  if(z < Hmax[x].x() && y < Hmax[x].y()) {
                    size_t gidx = offset(base.x() + x, y, z, imgSize.x(), imgSize.y());
                    if(gidx < imgMax)
                      data[gidx] = Hpix[idx];
                  }
                  idx++;
                }
          }
        }
      }
      catch(string& msg) {
        cerr << "voxelEditGPU::Error " << msg << endl;
        err = 1;
      }
      catch(std::exception& e) {
        cerr << "voxelEditGPU::Error" << e.what() << endl;
        err = 1;
      }
      catch(...) {
        cerr << "voxelEditGPU::Unknown exception" << endl;
        err = 1;
      }
    return (err);
  }
}
