/**============================================================================
 *
 * \file    run_artif_salt_pepper_noise.hpp
 * \brief   Process artificial data files with gaussian noise in them.
 *
 * \author  Leblanc Christophe
 * \date    04/04/2016
 * \e-mail: cleblancad@gmail.com
 *
 *===========================================================================*/

#ifndef TOMOPROCESS_RUN_ARTIF_SALT_PEPPER_NOISE_HPP
#define TOMOPROCESS_RUN_ARTIF_SALT_PEPPER_NOISE_HPP

#include "includes.hpp"

void run_artif_salt_pepper_noise()
{
  /*
   * A. Load images.
   */
  const unsigned int dimension(3); // 3D-image.
  std::ostream* outputMessageStream = &(std::cout); // Stream for messages.
  
  typedef typename itk::RGBPixel<unsigned char> InputPixelType;
  typedef float OutputPixelType;
  typedef Image<InputPixelType, OutputPixelType, dimension> ImageType;
  
  ImageType im(*outputMessageStream);
  std::string directory = "3DITK/artificial_salt_pepper_noise_50percent/";
  std::string extension = ".png";
  std::size_t fileStep(1), maxFilesToRead(100000);
  im.Read(directory, extension, fileStep, maxFilesToRead);
  im.SetVerboseOn();
  
  /*
   * B. Filters.
   */
  // 1. Streaming filter.
  /// Parameters.
  /// \numberOfStreamDivisions sets the number of pieces to divide the input. The upstream pipeline will be executed this many times.
  /// \createTemporaryFiles creates temporary files for each streaming slice.
  /// \removeTemporaryFiles removes the created temporary files.
  /// \os output stream for error messages.
  typedef StreamingFilter<ImageType::NextSelf> StreamingFilterType;
  typedef boost::shared_ptr<StreamingFilterType> StreamingFilterPtrType;
  
  StreamingFilterPtrType streamingFilter( new StreamingFilterType );
  const unsigned int numberOfStreamDivisions(2);
  const bool createTemporaryFiles = true, removeTemporaryFiles = true;
  
  streamingFilter->SetParameters(numberOfStreamDivisions, 
				 createTemporaryFiles, 
				 removeTemporaryFiles, 
				 *outputMessageStream);
  streamingFilter->ReleaseDataFlagOn();
  
  const StreamingFilterType::ImageInformationsContainerType* imageInformations = 
    streamingFilter->GetImageInformations();
  
  // 2. Spacing filter.
  /// Parameters.
  /// \spacing: new image spacing.
  typedef ChangeSpacingFilter<ImageType::NextSelf> ChangeSpacingFilterType;
  typedef boost::shared_ptr<ChangeSpacingFilterType> ChangeSpacingFilterPtrType;
  typedef ChangeSpacingFilterType::SpacingType SpacingType;
  
  ChangeSpacingFilterPtrType changeSpacingFilter( new ChangeSpacingFilterType );
  SpacingType spacing;
  
  // PHTR series (spacing in µm).
  spacing[0] = 1.0;
  spacing[1] = 1.0;
  spacing[2] = 1.0*fileStep; // 5.0*fileStep (for elongated version);
  
  changeSpacingFilter->SetDescription("changeSpacingFilter");
  changeSpacingFilter->SetParameters(spacing);
  changeSpacingFilter->ReleaseDataFlagOn();
  
  // 3. Automatic binary threshold filter.
  ///
  /// Parameters.
  /// \insideValue value for the foreground.
  /// \outsideValue value for the background.
  /// \numberOfHistogramBins number of histogram bins.
  /// \histogramBinMinimum minimum histogram bin value.
  /// \histogramBinMaximum maximum histogram bin value.
  typedef AutoBinaryThresholdFilter<ImageType::NextSelf> AutoBinaryThresholdFilterType;
  typedef boost::shared_ptr<AutoBinaryThresholdFilterType> AutoBinaryThresholdFilterPtrType;
  
  AutoBinaryThresholdFilterPtrType autoBinaryThresholdFilter( new AutoBinaryThresholdFilterType );
  const unsigned int numberOfHistogramBins(256);
  const AutoBinaryThresholdFilterType::InputPixelType histogramBinMinimum(0), 
    histogramBinMaximum(255);
  const AutoBinaryThresholdFilterType::OutputPixelType insideValue(0), outsideValue(1); // Image inverted, re-inverse here
  
  autoBinaryThresholdFilter->SetDescription("autoBinaryThresholdFilter");
  autoBinaryThresholdFilter->SetParameters(insideValue, outsideValue,
					   numberOfHistogramBins,
					   histogramBinMinimum, histogramBinMaximum);
  autoBinaryThresholdFilter->ReleaseDataFlagOn();
  
  
  // 4. Percolation filter.
  ///
  /// Parameters.
  /// \lookFor pixels' value to look for.
  /// \radius radii for the minimal cluster size (clusters of pixels smaller than that are deleted).
  /// \foregroundValue value for the selected pixels.
  /// \backgroundValue value for the deleted pixels.
  /// \useSpacing taking into account the image spacing for the computation of cluster sizes.
  typedef PercolationFilter<ImageType::NextSelf> PercolationFilterType;
  typedef boost::shared_ptr<PercolationFilterType> PercolationFilterPtrType;

  PercolationFilterPtrType percolationFilter( new PercolationFilterType );
  const bool percolationUseImageSpacing(true);
  PercolationFilterType::RadiusType percolationRadius;
  
  // Set radii in µm.
  percolationRadius[0] = 30;
  percolationRadius[1] = 30;
  percolationRadius[2] = 30;
  
  percolationFilter->SetDescription("percolationfilter");
  percolationFilter->SetParameters(outsideValue, percolationRadius, outsideValue, insideValue, 
				   percolationUseImageSpacing);
  percolationFilter->ReleaseDataFlagOn();
  
  //---------------------------------------------------
  
  typedef ErosionFilter<ImageType::NextSelf> ErosionFilterType;
  typedef boost::shared_ptr<ErosionFilterType> ErosionFilterPtrType;
  
  ErosionFilterPtrType erosionFilter( new ErosionFilterType );
  ErosionFilterType::RadiusType erosionRadius(1);
  
  erosionFilter->SetDescription("erosionFilter");
  erosionFilter->SetParameters(insideValue, outsideValue, erosionRadius);
  erosionFilter->ReleaseDataFlagOn();
  
  
  //---------------------------------------------------
  
  
  // 5. Distance filter.
  ///
  /// Parameters.
  /// \useSquaredDistance uses squared distance instead of euclidean distance.
  /// \useImageSpacing uses the image spacing for computing the distances.
  /// \background background value which defines the object. Usually this value is = 0.
  typedef DistanceFilter<ImageType::NextSelf> DistanceFilterType;
  typedef boost::shared_ptr<DistanceFilterType> DistanceFilterPtrType;
  
  const bool useSquaredDistance = false;
  const bool distanceUseImageSpacing = false; // Bug in maurer_distance_map_filter when this is true...
  DistanceFilterPtrType distanceFilter( new DistanceFilterType );
  
  distanceFilter->SetDescription("distanceFilter");
  distanceFilter->SetParameters(useSquaredDistance, 
				distanceUseImageSpacing, insideValue);
  distanceFilter->ReleaseDataFlagOn();
  
  // 6. Distance streaming correction filter.
  ///
  /// Parameter.
  /// \useImageSpacing uses the image spacing for computing the distances.
  typedef DistanceCorrectionFilter<ImageType::NextSelf> DistanceCorrectionFilterType;
  typedef boost::shared_ptr<DistanceCorrectionFilterType> DistanceCorrectionFilterPtrType;
  
  DistanceCorrectionFilterPtrType distanceCorrectionFilter( new DistanceCorrectionFilterType );
  
  distanceCorrectionFilter->SetDescription("distanceCorrectionFilter");
  distanceCorrectionFilter->SetParameters(distanceUseImageSpacing);
  distanceCorrectionFilter->ReleaseDataFlagOn();
  
  // 6. Maxima filter.
  ///
  /// Parameters.
  /// \foreground sets the value in the output image to consider as "foreground".
  /// \background1 sets the value in the output image to consider as "background".
  /// \background2 sets the value in the output image to consider as cell walls.
  /// \radius sets the radius of the considered neighborhood.
  /// \sigma value of scale for the derivative operator.
  /// \eigenTol tolerance factor for maxima selections.
  /// \maximumError maximum error for gaussian kernel calculation (usually 0.02 or less).
  /// \plateauTol two neighbors pixels with intensity difference less than PlateauTol are considered to be part of the same plateau.
  /// \useImageSpacing sets if we use the image spacing into the computations.
  /// \draw tells the filter if it has to draw into the output image.
  typedef MaximaFilter<ImageType::NextSelf> MaximaFilterType;
  typedef boost::shared_ptr<MaximaFilterType> MaximaFilterPtrType;
  
  MaximaFilterType::PixelType maximaBackground1(0), maximaBackground2(1), maximaForeground(2);
  unsigned int maximaRadius(5);
  double maximaSigma(1.0), maximaMaximumError(0.02), maximaPlateauTol(1.0e-5);
  double maximaEigenTol(0.1);
  bool maximaUseImageSpacing(true), maximaDraw(false);
  
  MaximaFilterPtrType maximaFilter( new MaximaFilterType );
  
  maximaFilter->SetDescription("maximaFilter");
  maximaFilter->SetParameters(maximaForeground, maximaBackground1, maximaBackground2,
    maximaRadius, maximaSigma, maximaEigenTol, maximaMaximumError, maximaPlateauTol,
    maximaUseImageSpacing, maximaDraw);
  maximaFilter->ReleaseDataFlagOn();
  
  MaximaFilterType::MaximaInformationsContainerType* maximaInformations =
    maximaFilter->GetMaximaInformations();
  
  //-------------------------------------------------------
  
  // For display purpose only.
  typedef DilationFilter<ImageType::NextSelf> DilationFilterType;
  typedef boost::shared_ptr<DilationFilterType> DilationFilterPtrType;
  
  DilationFilterType::RadiusType dilRadius(5);
  DilationFilterPtrType dilationFilter( new DilationFilterType );
  
  dilationFilter->SetDescription("dilationFilter");
  dilationFilter->SetParameters(maximaForeground, maximaBackground1, dilRadius); // For maxima
  dilationFilter->ReleaseDataFlagOn();
  
  //-------------------------------------------------------
  
  // 7. Centroids filter.
  ///
  /// Parameters.
  /// \foreground foreground pixels.
  /// \background backgroudn pixels.
  /// \feature feature pixel.
  /// \theta angular step (in radian) for evaluating points on an ellipsoidal surface. (Longitude).
  /// \phi angular step (in radian) for evaluating points on an ellipsoidal surface. (Latitude).
  /// \minNeighbours minimum number of neighbours to cluster a maximum.
  /// \overlappingRatio minimal volume overlapping ratio between two ellipsoids (\f$\in [0, 1]\f$).
  /// \maxTrials maximum number of trial points for estimating the surface overlapping ratio between two ellipsoids.
  /// \tolerance relative tolerance in the intersecting volume estimation between two ellipsoids.
  /// \maxInfos informations on maxima.
  /// \imagesInfos informations on images.
  /// \discardBoundaries sets if the maxima on boundaries should be discarded.
  /// \draw flag for drawing (true) or not (false) to the output image.
  /// \postExpand sets if a post-expansion of the ellipsoid is done after clustering.
  /// \computeMode sets the compute mode (full, fast).
  /// \computeQuality flag for computing the cell-ellipsoid matching quality.
  /// \verbose sets the verbose mode on (true) or off (false).
  /// \os output stream for messages.
  typedef CentroidsFilter<ImageType::NextSelf, MaximaFilterType::MaximaInformationsContainerType,
    StreamingFilterType::ImageInformationsContainerType> CentroidsFilterType;
  typedef boost::shared_ptr<CentroidsFilterType> CentroidsFilterPtrType;
  typedef typename CentroidsFilterType::ComputeModeType CentroidsComputeModeType;
    
  bool centroidsDraw(true), centroidsDiscardBoundaries(false), centroidQualityCompute(false);
  bool centroidsPostExpand(false);
  unsigned int centroidsComputeMode(CentroidsComputeModeType::full);
  typename CentroidsFilterType::InputPixelType centroidsFeature;
  typename CentroidsFilterType::OutputPixelType centroidsForeground(1), centroidsBackground(0);
  double centroidsTheta(M_PI/32.0), centroidsPhi(M_PI/32), centroidsOverlappingRatio(0.0);
  double centroidsTolerance(1.0e-1);
  unsigned int centroidsMinNeighbours(1), centroidsMaxTrials(4096);
  CentroidsFilterPtrType centroidsFilter( new CentroidsFilterType );
 
  if(maximaDraw) centroidsFeature = maximaBackground2;
  else centroidsFeature = 0.0;
  
  centroidsFilter->SetParameters(centroidsForeground,
				 centroidsBackground,
				 centroidsFeature,
				 centroidsTheta,
				 centroidsPhi,
				 centroidsMinNeighbours,
				 centroidsOverlappingRatio,
				 centroidsMaxTrials,
				 centroidsTolerance,
				 maximaInformations, 
				 imageInformations,
				 centroidsDiscardBoundaries,
				 centroidsDraw,
				 centroidsPostExpand,
				 centroidsComputeMode,
				 centroidQualityCompute,
				 im.GetVerbose(),
				 *outputMessageStream);
  centroidsFilter->SetDescription("centroidsFilter");
  centroidsFilter->ReleaseDataFlagOn();
  
  //-------------------------------------------------------
  
  // 8. Discretization filter.
  ///
  /// Parameters.
  /// \ellipsoidToleranceSize size tolerance for ellipsoids. Ellipsoids smaller than this will be erased.
  /// \maxInfos informations on maxima.
  /// \decimationMode decimation mode for small edges and faces (absolute or relative).
  /// \smallEdgeTolerance edge lengths below this tolerance will be decimated.
  /// \smallfaceTolerance face area below this tolerance will be decimated.
  /// \wantedFoamDensity wanted foam density (\$\in [0, 1]\$).
  /// \materialDensity density of wall's material of the foam.
  /// \fileName output file name.
  /// \meshScale characteristic length for the mesh.
  /// \discretizeTolerance tolerance for the optimization method.
  /// \intersectionTolerance tolerance for intersecting meshes.
  /// \moveTolerance amount by which a point of an intersecting mesh will be moved.
  /// \maxIter maximum number of allowed iterations for the optimization method.
  /// \imagesInfos informations on images.
  /// \verbose verbose mode.
  /// \outputStream output stream for verbosity.
  typedef DiscretizeFilter<ImageType::NextSelf, MaximaFilterType::MaximaInformationsContainerType,
    StreamingFilterType::ImageInformationsContainerType> DiscretizeFilterType;
  typedef boost::shared_ptr<DiscretizeFilterType> DiscretizeFilterPtrType;
  
  int decimationMode = DiscretizeFilterType::absolute;
  double ellipsoidToleranceSize(0.0), smallEdgeTolerance(1.0e-2), smallfaceTolerance(1.0e-2),
    wantedFoamDensity(0.0), meshScale(50.0), materialDensity(0.0), discretizeTolerance(1.0e-3),
    intersectionTolerance(1.0e-5), moveTolerance(0.1);
  unsigned int discretizeMaxIter(100);
  std::string discretizeFileName("test.geo");
  
  DiscretizeFilterPtrType discretizeFilter( new DiscretizeFilterType );
  
  discretizeFilter->SetParameters(ellipsoidToleranceSize,
                                  maximaInformations,
				  decimationMode,
				  smallEdgeTolerance,
				  smallfaceTolerance,
				  wantedFoamDensity,
				  materialDensity,
				  discretizeFileName,
				  meshScale,
				  discretizeTolerance,
				  intersectionTolerance,
				  moveTolerance,
				  discretizeMaxIter,
				  imageInformations,
				  im.GetVerbose(),
				  *outputMessageStream);
  discretizeFilter->SetDescription("discretizeFilter");
  discretizeFilter->ReleaseDataFlagOn();
  
  /*
   * C. Output.
   */
  try
  {
    // Changes image spacing
    // (the previous filters do not use spacing information).
    im << changeSpacingFilter;
    
    // Applies a binray threshold.
    im << autoBinaryThresholdFilter;
    
    // Applies a cleaning.
    //im << erosionFilter;
    im << percolationFilter;
    
    // Applies a distance transform.
    im << distanceFilter;
    
    // Computes the streaming correction to the distance map.
    im << distanceCorrectionFilter;
    
    // Locates local maxima.
    im << maximaFilter;
    //im << dilationFilter;
    
    // Streaming filter.
    //im << streamingFilter;
    
    // Centroids filter.
    im << centroidsFilter;
    
    // Discretization filter.
    //im << discretizeFilter;
    
    // Write file.
    im.Write("artif_salt_pepper_noise_50percent/F_centroids_full1.vtk");
    //im.Update();
  }
  catch(itk::ExceptionObject &err)
  {
    (*outputMessageStream) << err << std::endl;
  }
}

#endif // TOMOPROCESS_RUN_ARTIF_SALT_PEPPER_NOISE_HPP