//
// 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 "AutoTrim.hpp"
#include "Information.hpp"
#include "Progress.hpp"

namespace mgx 
{
  bool AutoTrim::run(Stack* stack, Store* store, int threshold)
  {
    Point3u size = stack->size();
    Point3u pmin(stack->size().x() + stack->size().y() + stack->size().z()), pmax(0u);
    progressStart("Auto-trim stack", 2 * size.y() * size.z());
    int adv = 0, adv_step = 2 * size.y() * size.z() / 100;
    const HVecUS& data = store->data();
    SETSTATUS("Search bounds.");
    for(uint z = 0, k = 0; z < size.z(); ++z)
      for(uint y = 0; y < size.y(); ++y, ++adv) {
        for(uint x = 0; x < size.x(); ++x, ++k) {
          if(data[k] > threshold) {
            if(x < pmin.x())
              pmin.x() = x;
            if(y < pmin.y())
              pmin.y() = y;
            if(z < pmin.z())
              pmin.z() = z;
            if(x > pmax.x())
              pmax.x() = x;
            if(y > pmax.y())
              pmax.y() = y;
            if(z > pmax.z())
              pmax.z() = z;
          }
        }
        if((adv % adv_step == 0)and !progressAdvance(adv))
          userCancel();
      }
    Point3u newSize = (pmax - pmin) + 1u;
    if(newSize == size) {
      SETSTATUS("No change done, the volume already occupies the whole space.");
      return true;
    }
    SETSTATUS("Copy data.");
    unsigned int totSize = newSize.x() * newSize.y() * newSize.z();
    if(totSize == 0)
      throw QString("Error, found a volume of 0, is the image empty?");
    HVecUS mainData(totSize);
    HVecUS workData(totSize);
    unsigned int stride_y = (size.y() - newSize.y()) * size.x();
    ushort* srcMain = &stack->main()->data()[stack->offset(Point3i(pmin))];
    ushort* dstMain = &mainData[0];
    ushort* srcWork = &stack->work()->data()[stack->offset(Point3i(pmin))];
    ushort* dstWork = &workData[0];
    for(uint z = 0; z < newSize.z(); ++z) {
      for(uint y = 0; y < newSize.y(); ++y, ++adv) {
        memcpy(dstMain, srcMain, sizeof(ushort) * newSize.x());
        dstMain += newSize.x();
        srcMain += size.x();
  
        memcpy(dstWork, srcWork, sizeof(ushort) * newSize.x());
        dstWork += newSize.x();
        srcWork += size.x();
      }
      srcMain += stride_y;
      srcWork += stride_y;
      if((adv % adv_step == 0)and !progressAdvance(adv))
        userCancel();
    }
  
    if(!progressAdvance(2 * size.y() * size.z()))
      userCancel();
    SETSTATUS("Resize actual stack.");
  
    // Compute new origin
    stack->setSize(newSize);
  
    stack->setOrigin(stack->imageToWorld(Point3f(pmin) - Point3f(.5, .5, .5)));
  
    stack->main()->data() = mainData;
    stack->work()->data() = workData;
  
    stack->main()->changed();
    stack->work()->changed();
    return true;
  }
  
  REGISTER_PROCESS(AutoTrim);
}
