//
// 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.
// 
#ifndef STACKPROCESSFILTERS_HPP
#define STACKPROCESSFILTERS_HPP

#include <Process.hpp>

namespace mgx 
{
  ///\addtogroup Process
  ///@{
  /**
   * \class AverageStack StackProcessFilters.hpp <StackProcessFilters.hpp>
   *
   * Compute the local average of each voxel of the stack.
   *
   * More precisely, this module perform the convolution of the stack with the
   * square kernel of given radius.
   */
  class mgxBase_EXPORT AverageStack : public Process 
  {
  public:
    AverageStack(const Process& process) : Process(process) 
    {
      setName("Stack/Filters/Average");
      setDesc("Average stack data with a square filter.");
      setIcon(QIcon(":/images/AvgPix.png"));
      
      addParm("X Radius","X Radius","1");
      addParm("Y Radius","Y Radius","1");
      addParm("Z Radius","Z Radius","1");
      addParm("Steps","Steps","1");						
	}
  
    bool run()
    {
      if(!checkState().store(STORE_NON_EMPTY))
        return false;
      Store* input = currentStack()->currentStore();
      Store* output = currentStack()->work();
      Point3i r(parm("X Radius").toInt(), parm("Y Radius").toInt(), parm("Z Radius").toInt());
      bool res = run(input, output, r, parm("Steps").toUInt());
      if(res) {
        input->hide();
        output->show();
      }
      return res;
    }
  
    bool run(const Store* input, Store* output, Point3i radius, uint steps);
  };
  
  /**
   * Uses quickselect to compute the local median of each voxel of the stack.
   */
  class mgxBase_EXPORT MedianStack : public Process 
  {
  public:
    MedianStack(const Process& process) : Process(process) 
    {
      setName("Stack/Filters/Median");
      setDesc("Median of stack data.");
      setIcon(QIcon(":/images/AvgPix.png"));
      
      addParm("X Radius","X Radius","1");
      addParm("Y Radius","Y Radius","1");
      addParm("Z Radius","Z Radius","1");				
	}
  
    bool run()
    {
      if(!checkState().store(STORE_NON_EMPTY))
        return false;
      Store* input = currentStack()->currentStore();
      Store* output = currentStack()->work();
      Point3i r(parm("X Radius").toInt(), parm("Y Radius").toInt(), parm("Z Radius").toInt());
      bool res = run(input, output, r);
      if(res) {
        input->hide();
        output->show();
      }
      return res;
    }
  
    bool run(const Store* input, Store* output, Point3i radius); 
  };
  
  /**
   * Sepearable median-of-medians approximation to median filter
   */
  class mgxBase_EXPORT MedianOfMediansStack : public Process 
  {
  public:
    MedianOfMediansStack(const Process& process) : Process(process) 
    {
		
      setName("Stack/Filters/Median of Medians");
      setDesc("Median of medians of stack data.");
      setIcon(QIcon(":/images/AvgPix.png"));
      
      addParm("X Radius","X Radius","1");
      addParm("Y Radius","Y Radius","1");
      addParm("Z Radius","Z Radius","1");				
	}  
	
    bool run()
    {
      if(!checkState().store(STORE_NON_EMPTY))
        return false;
      Store* input = currentStack()->currentStore();
      Store* output = currentStack()->work();
      Point3i r(parm("X Radius").toInt(), parm("Y Radius").toInt(), parm("Z Radius").toInt());
      bool res = run(input, output, r);
      if(res) {
        input->hide();
        output->show();
      }
      return res;
    }
 
    bool run(const Store* input, Store* output, Point3i radius);

  };
  
  /**
   * \class BrightenStack StackProcessFilters.hpp <StackProcessFilters.hpp>
   *
   * Change the global luminosity of the stack. The value of each voxel will be
   * multiplied by the brightening constant. So a value > 1 will brighten the
   * stack, while a value < 1 will darken the stack
   *
   * Useful examples:
   *  - to convert a stack from 12bits to 16bits, use 16;
   *  - to convert a stack from 8bits to 16bits, use 256.
   */
  class mgxBase_EXPORT BrightenStack : public Process {
  public:
    BrightenStack(const Process& process) : Process(process) 
    {
      setName("Stack/Filters/Brighten Darken");
      setDesc("Brighten or Darken stack.\n"
      "A value > 1 will brighten the stack, < 1 will darken it.\n"
      "To convert from a 12bit to a 16bit stack, use 16.\n"
      "To convert from a 8bit to a 16bit stack, use 256.");
      setIcon(QIcon(":/images/Brightness.png"));
      
      addParm("Amount","Amount to multiply voxels. Use 16 for 12bit images and 256 for 8bit images","16.0");					
	}
  
    bool run()
    {
      if(!checkState().store(STORE_NON_EMPTY))
        return false;
      Stack* s = currentStack();
      Store* input = s->currentStore();
      Store* output = s->work();
      bool res = run(input, output, parm("Amount").toFloat());
      if(res) {
        input->hide();
        output->show();
      }
      return res;
    }
    
    bool run(const Store* input, Store* output, float brightness);
  
  };
  
  /**
   * \class BinarizeStack StackProcessFilters.hpp <StackProcessFilters.hpp>
   *
   * Binarize the stack: set the intensity of any voxel greater than the
   * threshold to 65535 and all other voxels to 0.
   */
  class mgxBase_EXPORT BinarizeStack : public Process 
  {
  public:
    BinarizeStack(const Process& process) : Process(process) 
    {
      setName("Stack/Filters/Binarize");
      setDesc("Transform the stack to binary (65535 or 0).\n"
      "All voxels with an intensity greater than the threshold to 65535, while to others are set to 0.");
      setIcon(QIcon(":/images/Brightness.png"));
      
      addParm("Threshold","Voxels above this threshold will be assigned 1.","5000");				
	}
  
    bool run()
    {
      if(!checkState().store(STORE_NON_EMPTY))
        return false;
      Stack* s = currentStack();
      Store* input = s->currentStore();
      Store* output = s->work();
      bool res = run(input, output, parm("Threshold").toUShort());
      if(res) {
        input->hide();
        output->show();
      }
      return res;
    }
    
    bool run(const Store* input, Store* output, ushort threshold);
  
  };
  
  /**
   * \class ColorGradient StackProcessFilters.hpp <StackProcessFilters.hpp>
   *
   * Compute the gradient of the image in the Z axis of the image.
   */
  class mgxBase_EXPORT ColorGradient : public Process 
  {
  public:
    ColorGradient(const Process& process) : Process(process) 
    {
      setName("Stack/Filters/Color Gradient");
      setDesc("Compute color gradient in Z direction");
      setIcon(QIcon(":/images/ColorGrad.png"));
      
      addParm("Z Divisor","Factor by which the gradient is divided by.","5.0");				
	}
  
    bool run()
    {
      if(!checkState().store(STORE_NON_EMPTY))
        return false;
      Stack* s = currentStack();
      Store* input = s->currentStore();
      Store* output = s->work();
      bool res = run(input, output, parm("Z Divisor").toFloat());
      if(res) {
        input->hide();
        output->show();
      }
      return res;
    }
    
    bool run(const Store* input, Store* output, float colorGradDiv);
  };
  
  /**
   * \class InvertStack StackProcessFilters.hpp <StackProcessFilters.hpp>
   *
   * Invert the intensity of all voxels in the stack.
   */
  class mgxBase_EXPORT InvertStack : public Process 
  {
  public:
    InvertStack(const Process& process) : Process(process) 
    {
      setName("Stack/Filters/Invert");
      setDesc("Invert the stack");	
      setIcon(QIcon(":/images/Invert.png"));			
	}
  
    bool run()
    {
      if(!checkState().store(STORE_NON_EMPTY))
        return false;
      Stack* s = currentStack();
      Store* input = s->currentStore();
      Store* output = s->work();
      bool res = run(input, output);
      if(res) {
        input->hide();
        output->show();
      }
      return res;
    }
    
    bool run(const Store* input, Store* output);
  
  };
  
  /**
   * \class FilterStack StackProcessFilters.hpp <StackProcessFilters.hpp>
   *
   * Clip the intensity of the stack to the interval [Low Threshold, High Threshold].
   */
  class mgxBase_EXPORT FilterStack : public Process 
  {
  public:
    FilterStack(const Process& process) : Process(process) 
    {
      setName("Stack/Filters/Trim High Low Values");
      setDesc("Clip the voxel intensities to the interval [Low Threshold, High Threshold].");
      setIcon(QIcon(":/images/Filter.png"));
      
      addParm("Low Threshold","Lower bound","1000");
      addParm("High Threshold","Upper bound","0");				
	}
  
    bool run()
    {
      if(!checkState().store(STORE_NON_EMPTY))
        return false;
      Stack* s = currentStack();
      Store* input = s->currentStore();
      Store* output = s->work();
      bool res = run(input, output, parm("Low Threshold").toUInt(), parm("High Threshold").toUInt());
      if(res) {
        input->hide();
        output->show();
      }
      return res;
    }
    
    bool run(const Store* input, Store* output, uint lowFilter, uint highFilter);
  
  };
  
  /**
   * \class GaussianBlurStack StackProcessFilters.hpp <StackProcessFilters.hpp>
   *
   * Blur the stack using a Gaussian kernel or specified radius.
   *
   * The radius given correspond to the standard deviation of the Gaussian
   * kernel along the various axes. The kernel itself will be 3 time as big.
   */
  class mgxBase_EXPORT GaussianBlurStack : public Process 
  {
  public:
    GaussianBlurStack(const Process& process) : Process(process) 
    {
      setName("Stack/Filters/Gaussian Blur Stack");
      setDesc("Blur the stack, radius = 3 x Sigma");
      setIcon(QIcon(":/images/Gaussian.png"));
      
      addParm("X Sigma (µm)","X Sigma (µm)","0.3");
      addParm("Y Sigma (µm)","Y Sigma (µm)","0.3");
      addParm("Z Sigma (µm)","Z Sigma (µm)","0.3");				
	}
  
    bool run()
    {
      if(!checkState().store(STORE_NON_EMPTY))
        return false;
      Stack* s = currentStack();
      Store* input = s->currentStore();
      Store* output = s->work();
      Point3f sigma(parm("X Sigma (µm)").toFloat(), parm("Y Sigma (µm)").toFloat(), parm("Z Sigma (µm)").toFloat());
      bool res = run(input, output, sigma);
      if(res) {
        input->hide();
        output->show();
      }
      return res;
    }
    
    bool run(const Store* input, Store* output, Point3f sigma);
  };

  /**
   * \class DiffGaussiansStack StackProcessFilters.hpp
   *
   * Calculate the difference of Gaussians, i.e. first - second.
   *
   * The radius given correspond to the standard deviation of the Gaussian
   * kernels along the various axes. The kernel itself will be 3 times as big.
   */
  class mgxBase_EXPORT DiffGaussians : public Process 
  {
  public:
    DiffGaussians(const Process& process) : Process(process) 
    {
      setName("Stack/Filters/Diff Gaussians");
      setDesc("Difference of Gaussians, first - second");
      setIcon(QIcon(":/images/DiffGaussians.png"));
      
      addParm("X1 Sigma (µm)","X1 Sigma (µm)","0.3");
      addParm("Y1 Sigma (µm)","Y1 Sigma (µm)","0.3");
      addParm("Z1 Sigma (µm)","Z1 Sigma (µm)","0.3");
      addParm("X2 Sigma (µm)","X2 Sigma (µm)","10.0");
      addParm("Y2 Sigma (µm)","Y2 Sigma (µm)","10.0");
      addParm("Z2 Sigma (µm)","Z2 Sigma (µm)","10.0");				
	}  
	
    bool run()
    {
      if(!checkState().store(STORE_NON_EMPTY))
        return false;
      Stack* s = currentStack();
      Store* input = s->currentStore();
      Store* output = s->work();
      Point3f sigma1(parm("X1 Sigma (µm)").toFloat(), parm("Y1 Sigma (µm)").toFloat(), parm("Z1 Sigma (µm)").toFloat());
      Point3f sigma2(parm("X2 Sigma (µm)").toFloat(), parm("Y2 Sigma (µm)").toFloat(), parm("Z2 Sigma (µm)").toFloat());
      bool res = run(input, output, sigma1, sigma2);
      if(res) {
        input->hide();
        output->show();
      }
      return res;
    }
    
    bool run(const Store* input, Store* output, Point3f sigma1, Point3f sigma2);

  };

  /**
   * \class SharpenStack StackProcessFilters.hpp <StackProcessFilters.hpp>
   *
   * Apply an "unsharp" filter to the stack.
   *
   * The unsharp filter consists in subtracting from the image the blured
   * version of itself. The parameters of the filter correspond to the blurring
   * and a multiplicating factor for the blured image.
   */
  class mgxBase_EXPORT SharpenStack : public Process 
  {
  public:
    SharpenStack(const Process& process) : Process(process) 
    {
      setName("Stack/Filters/Sharpen Stack");
      setDesc("Sharpen the stack, radius = 3 x Sigma");
      setIcon(QIcon(":/images/Sharpen.png"));
      
      addParm("X Sigma (µm)","X Sigma (µm)","0.2");
      addParm("Y Sigma (µm)","Y Sigma (µm)","0.2");
      addParm("Z Sigma (µm)","Z Sigma (µm)","1.0");
      addParm("Amount","Amount","1.0");
      addParm("Passes","Passes","1");				
	}
  
    bool run()
    {
      if(!checkState().store(STORE_NON_EMPTY))
        return false;
      Stack* s = currentStack();
      Store* input = s->currentStore();
      Store* output = s->work();
      Point3f sigma(parm("X Sigma (µm)").toFloat(), parm("Y Sigma (µm)").toFloat(), parm("Z Sigma (µm)").toFloat());
      bool res = run(input, output, sigma, parm("Amount").toFloat(), parm("Passes").toInt());
      if(res) {
        input->hide();
        output->show();
      }
      return res;
    }
    
    bool run(const Store* input, Store* output, 
                                           const Point3f &sigma, float amount, int passes);
  };
  
  /**
   * \class ApplyKernelStack StackProcessFilters.hpp <StackProcessFilters.hpp>
   *
   * Apply a user-defined convolution kernel to the current image.
   *
   * The kernel must be separable, so only 1D kernels for the X, Y and Z axes
   * are specified.
   */
  //class mgxBase_EXPORT ApplyKernelStack : public Process {
  //public:
    //ApplyKernelStack(const Process& process) : Process(process) 
    //{
      //setName("Stack/Filters/Apply Separable Kernel");
      //setDesc("Kernel must have odd number of values");
      //setIcon(QIcon(":/images/Kernel.png"));
      
      //addParm("X Kernel","X Kernel","-.1 1.2 -.1");
      //addParm("Y Kernel","Y Kernel","-.1 1.2 -.1");
      //addParm("Z Kernel","Z Kernel","-.1 1.2 -.1");	
	//}
  
    //bool run()
    //{
      //if(!checkState().store(STORE_NON_EMPTY))
        //return false;
      //Stack* s = currentStack();
      //Store* input = s->currentStore();
      //Store* output = s->work();
      //HVecF kernels[3];
      //bool ok;
      //QString fields[3] = { "x", "y", "z" };
      //for(int i = 0; i < 3; i++) {
        //QString str = QString(parms[i]).trimmed();
        //QStringList field_values = str.split(" ", QString::SkipEmptyParts);
        //kernels[i].resize(field_values.size());
        //for(size_t j = 0; j < (size_t)field_values.size(); ++j) {
          //kernels[i][j] = field_values[j].toFloat(&ok);
          //if(not ok)
            //return setErrorMessage(
              //QString("Value %1 of field %2 if not a valid number.").arg(j).arg(fields[i]));
        //}
      //}
  
      //bool res = run(input, output, kernels[0], kernels[1], kernels[2]);
      //if(res) {
        //input->hide();
        //output->show();
      //}
      //return res;
    //}
    //bool run(const Store* input, Store* output, const HVecF& kernelX, const HVecF& kernelY,
                    //const HVecF& kernelZ);
  //};
  
  /**
   * \class NormalizeStack StackProcessFilters.hpp <StackProcessFilters.hpp>
   *
   * Normalize the intensity of the stack.
   *
   * This process first dilates, then blur the image to obtain a local maximum
   * voxel intensity. The value of the current voxel is then divided by this
   * local maximum.
   */
  class mgxBase_EXPORT NormalizeStack : public Process 
  {
  public:
    NormalizeStack(const Process& process) : Process(process) 
    {
      setName("Stack/Filters/Normalize Stack");
      setDesc("Normalize the stack");
      setIcon(QIcon(":/images/Normalize.png"));
      
      addParm("X Radius (µm)","X Radius (µm) for the locality of normalization","5.0");
      addParm("Y Radius (µm)","Y Radius (µm) for the locality of normalization","5.0");
      addParm("Z Radius (µm)","Z Radius (µm) for the locality of normalization","5.0");
      addParm("X Sigma (µm)","X Sigma (µm) for pre-blur","0.5");
      addParm("Y Sigma (µm)","Y Sigma (µm) for pre-blur","0.5");
      addParm("Z Sigma (µm)","Z Sigma (µm) for pre-blur","0.5");
      addParm("Threshold","Threshold under which pixels are cleared considered as background","10000");
      addParm("Blur factor","Relative contribution of blurred vs unblurred dilation","0.7");						
	}
  
    bool run()
    {
      if(!checkState().store(STORE_NON_EMPTY))
        return false;
      Stack* s = currentStack();
      Store* input = s->currentStore();
      Store* output = s->work();
      Point3f radius(parm("X Radius (µm)").toFloat(), parm("Y Radius (µm)").toFloat(), parm("Z Radius (µm)").toFloat());
      Point3f sigma(parm("X Sigma (µm)").toFloat(), parm("Y Sigma (µm)").toFloat(), parm("Z Sigma (µm)").toFloat());
      bool res = run(input, output, radius, sigma, parm("Threshold").toUInt(), parm("Blur factor").toFloat());
      if(res) {
        input->hide();
        output->show();
      }
      return res;
    }
    
    bool run(const Store* input, Store* output, Point3f radius, Point3f sigma, 
           uint threshold, float blurFactor);
  };
  
  /**
   * \class CImgGaussianBlurStack StackProcessFilters.hpp <StackProcessFilters.hpp>
   *
   * CImg implementation of the gaussian blur
   */
  class mgxBase_EXPORT CImgGaussianBlurStack : public Process 
  {
  public:
    CImgGaussianBlurStack(const Process& process) : Process(process) 
    {
      setName("Stack/CImage/Gaussian Blur");
      setDesc("CImage Gaussian Blur");
      setIcon(QIcon(":/images/Blur.png"));
      
      addParm("Radius","Radius","5");						
	}
  
    bool run()
    {
      if(!checkState().store(STORE_NON_EMPTY))
        return false;
      Stack* s = currentStack();
      Store* input = s->currentStore();
      Store* output = s->work();
      bool res = run(input, output, parm("Radius").toUInt());
      if(res) {
        input->hide();
        output->show();
      }
      return res;
    }
  
    bool run(const Store* input, Store* output, uint radius);
  };
	
    /**
   * \class CImgMedianBlurStack StackProcessFilters.hpp <StackProcessFilters.hpp>
   *
   * CImg implementation of the median blur
   */
  class mgxBase_EXPORT CImgMedianBlurStack : public Process 
  {
  public:
    CImgMedianBlurStack(const Process& process) : Process(process) 
    {
      setName("Stack/CImage/Median Blur");
      setDesc("CImage Median Blur");
      setIcon(QIcon(":/images/Blur.png"));
      
      addParm("Radius (voxels)","Size of the median filter.","3");				
	}
  
    bool run()
    {
      if(!checkState().store(STORE_NON_EMPTY))
        return false;
      Stack* s = currentStack();
      Store* input = s->currentStore();
      Store* output = s->work();
      bool res = run(input, output, parm("Radius (voxels)").toUInt());
      if(res) {
        input->hide();
        output->show();
      }
      return res;
    }
  
    bool run(const Store* input, Store* output, uint radius);
  };

  /**
   * \class CImgLaplaceStack StackProcessFilters.hpp <StackProcessFilters.hpp>
   *
   * Compute the Laplacian of the stack (using CImg).
   */
  class mgxBase_EXPORT CImgLaplaceStack : public Process 
  {
  public:
    CImgLaplaceStack(const Process& process) : Process(process) 
    {
      setName("Stack/CImage/Laplace Transform"); 
      setIcon(QIcon(":/images/Laplace.png"));		
	}
  
    bool run()
    {
      if(!checkState().store(STORE_NON_EMPTY))
        return false;
      Stack* s = currentStack();
      Store* input = s->currentStore();
      Store* output = s->work();
      bool res = run(input, output);
      if(res) {
        input->hide();
        output->show();
      }
      return res;
    }
    
    bool run(const Store* input, Store* output);
  };

   /**
   * \class TopHatStack StackProcessFilters.hpp <StackProcessFilters.hpp>
   *
   *
   * Perform the top-hat transformation on a stack. This filter is often used
   * for background subtraction
   */
  class mgxBase_EXPORT TopHatStack : public Process 
  {
  public:
    TopHatStack(const Process& process) : Process(process) 
    {
      setName("Stack/Filters/Top Hat");
      setDesc("Perform Top Hat transformation on stack, used for background subtraction");
      setIcon(QIcon(":/images/TopHat.png"));
      
      addParm("X Radius","X Radius","10");
      addParm("Y Radius","Y Radius","10");
      addParm("Z Radius","Z Radius","10");
      addParm("Round Nhbd","Use a round (Ellipsoidal) neighborhood","No",booleanChoice());					
	}
  
    bool run()
    {
      if(!checkState().store(STORE_NON_EMPTY))
        return false;
      Stack* s = currentStack();
      Store* input = s->currentStore();
      Store* output = s->work();
      bool res = run(s, input, output, 
              parm("X Radius").toUInt(), parm("Y Radius").toUInt(), parm("Z Radius").toUInt(), stringToBool(parm("Round Nhbd")));
      if(res) {
        input->hide();
        output->show();
      }
      return res;
    }
    
    bool run(Stack* stack, const Store* input, Store* output, 
                   uint xradius, uint yradius, uint zradius, bool roundNhbd); 
  }; 

   /**
   * \class AddNoiseStack StackProcessFilters.hpp <StackProcessFilters.hpp>
   *
   * Add noise to the stack
   */
  class mgxBase_EXPORT AddNoiseStack : public Process 
  {
  public:
    AddNoiseStack(const Process& process) : Process(process) 
    {
      setName("Stack/Filters/Add Noise");
      setDesc("Perform Top Hat transformation on stack, used for background subtraction");
      setIcon(QIcon(":/images/TopHat.png"));
      
      addParm("Amount","Amount of noise to add","1000");
	  }
  
    bool run()
    {
      if(!checkState().store(STORE_NON_EMPTY))
        return false;
      Stack* s = currentStack();
      Store* input = s->currentStore();
      Store* output = s->work();
      bool res = run(s, input, output, parm("Amount").toUInt());
      if(res) {
        input->hide();
        output->show();
      }
      return res;
    }
    
    bool run(Stack* stack, const Store* input, Store* output, uint amount); 
  }; 

  ///@}
}

#endif 
