/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         https://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/
#ifndef itkStatisticsLabelObject_h
#define itkStatisticsLabelObject_h

#include "itkShapeLabelObject.h"
#include "itkHistogram.h"

namespace itk
{
/**
 * \class StatisticsLabelObject
 *  \brief A Label object to store the common attributes related to the statistics of the object
 *
 * StatisticsLabelObject stores  the common attributes related to the statistics of the object
 *
 * \author Gaetan Lehmann. Biologie du Developpement et de la Reproduction, INRA de Jouy-en-Josas, France.
 *
 * This implementation was taken from the Insight Journal paper:
 * https://www.insight-journal.org/browse/publication/176
 *
 * \ingroup DataRepresentation
 * \ingroup ITKLabelMap
 */
template <typename TLabel, unsigned int VImageDimension>
class ITK_TEMPLATE_EXPORT StatisticsLabelObject : public ShapeLabelObject<TLabel, VImageDimension>
{
public:
  ITK_DISALLOW_COPY_AND_MOVE(StatisticsLabelObject);

  /** Standard class type aliases */
  using Self = StatisticsLabelObject;
  using Superclass = ShapeLabelObject<TLabel, VImageDimension>;
  using LabelObjectType = typename Superclass::LabelObjectType;
  using Pointer = SmartPointer<Self>;
  using ConstPointer = SmartPointer<const Self>;
  using ConstWeakPointer = WeakPointer<const Self>;

  /** Method for creation through the object factory. */
  itkNewMacro(Self);

  /** \see LightObject::GetNameOfClass() */
  itkOverrideGetNameOfClassMacro(StatisticsLabelObject);

  using LabelMapType = LabelMap<Self>;

  static constexpr unsigned int ImageDimension = VImageDimension;

  using typename Superclass::IndexType;

  using PointType = Point<double, Self::ImageDimension>;

  using LabelType = TLabel;

  using typename Superclass::LineType;

  using typename Superclass::LengthType;

  using MatrixType = Matrix<double, Self::ImageDimension, Self::ImageDimension>;

  using VectorType = Vector<double, Self::ImageDimension>;

  using HistogramType = Statistics::Histogram<double>;

  using typename Superclass::AttributeType;
  static constexpr AttributeType MINIMUM = 200;
  static constexpr AttributeType MAXIMUM = 201;
  static constexpr AttributeType MEAN = 202;
  static constexpr AttributeType SUM = 203;
  static constexpr AttributeType STANDARD_DEVIATION = 204;
  static constexpr AttributeType VARIANCE = 205;
  static constexpr AttributeType MEDIAN = 206;
  static constexpr AttributeType MAXIMUM_INDEX = 207;
  static constexpr AttributeType MINIMUM_INDEX = 208;
  static constexpr AttributeType CENTER_OF_GRAVITY = 209;
  //  static constexpr AttributeType CENTRAL_MOMENTS = 210;
  static constexpr AttributeType WEIGHTED_PRINCIPAL_MOMENTS = 211;
  static constexpr AttributeType WEIGHTED_PRINCIPAL_AXES = 212;
  static constexpr AttributeType KURTOSIS = 213;
  static constexpr AttributeType SKEWNESS = 214;
  static constexpr AttributeType WEIGHTED_ELONGATION = 215;
  static constexpr AttributeType HISTOGRAM = 216;
  static constexpr AttributeType WEIGHTED_FLATNESS = 217;

  static AttributeType
  GetAttributeFromName(const std::string & s)
  {
    if (s == "Minimum")
    {
      return MINIMUM;
    }
    else if (s == "Maximum")
    {
      return MAXIMUM;
    }
    else if (s == "Mean")
    {
      return MEAN;
    }
    else if (s == "Sum")
    {
      return SUM;
    }
    else if (s == "StandardDeviation")
    {
      return STANDARD_DEVIATION;
    }
    else if (s == "Variance")
    {
      return VARIANCE;
    }
    else if (s == "Median")
    {
      return MEDIAN;
    }
    else if (s == "MaximumIndex")
    {
      return MAXIMUM_INDEX;
    }
    else if (s == "MinimumIndex")
    {
      return MINIMUM_INDEX;
    }
    else if (s == "CenterOfGravity")
    {
      return CENTER_OF_GRAVITY;
    }
    /*
    else if( s == "CentralMoments" )
      {
      return CENTRAL_MOMENTS;
      }
    */
    else if (s == "WeightedPrincipalMoments")
    {
      return WEIGHTED_PRINCIPAL_MOMENTS;
    }
    else if (s == "WeightedPrincipalAxes")
    {
      return WEIGHTED_PRINCIPAL_AXES;
    }
    else if (s == "Kurtosis")
    {
      return KURTOSIS;
    }
    else if (s == "Skewness")
    {
      return SKEWNESS;
    }
    else if (s == "WeightedElongation")
    {
      return WEIGHTED_ELONGATION;
    }
    else if (s == "Histogram")
    {
      return HISTOGRAM;
    }
    else if (s == "WeightedFlatness")
    {
      return WEIGHTED_FLATNESS;
    }
    // can't recognize the name
    return Superclass::GetAttributeFromName(s);
  }

  static std::string
  GetNameFromAttribute(const AttributeType & a)
  {
    switch (a)
    {
      case MINIMUM:
        return "Minimum";
      case MAXIMUM:
        return "Maximum";
      case MEAN:
        return "Mean";
      case SUM:
        return "Sum";
      case STANDARD_DEVIATION:
        return "StandardDeviation";
      case VARIANCE:
        return "Variance";
      case MEDIAN:
        return "Median";
      case MAXIMUM_INDEX:
        return "MaximumIndex";
      case MINIMUM_INDEX:
        return "MinimumIndex";
      case CENTER_OF_GRAVITY:
        return "CenterOfGravity";
      /*      case CENTRAL_MOMENTS:
              return "CentralMoments";*/
      case WEIGHTED_PRINCIPAL_MOMENTS:
        return "WeightedPrincipalMoments";
      case WEIGHTED_PRINCIPAL_AXES:
        return "WeightedPrincipalAxes";
      case KURTOSIS:
        return "Kurtosis";
      case SKEWNESS:
        return "Skewness";
      case WEIGHTED_ELONGATION:
        return "WeightedElongation";
      case HISTOGRAM:
        return "Histogram";
      case WEIGHTED_FLATNESS:
        return "WeightedFlatness";
    }
    // can't recognize the name
    return Superclass::GetNameFromAttribute(a);
  }

  using RegionType = ImageRegion<Self::ImageDimension>;

  using typename Superclass::CentroidType;

  template <typename TSourceLabelObject>
  void
  CopyAttributesFrom(const TSourceLabelObject * src)
  {
    Superclass::template CopyAttributesFrom<TSourceLabelObject>(src);

    m_Minimum = src->GetMinimum();
    m_Maximum = src->GetMaximum();
    m_Mean = src->GetMean();
    m_Sum = src->GetSum();
    m_StandardDeviation = src->GetStandardDeviation();
    m_Variance = src->GetVariance();
    m_Median = src->GetMedian();
    m_MaximumIndex = src->GetMaximumIndex();
    m_MinimumIndex = src->GetMinimumIndex();
    m_CenterOfGravity = src->GetCenterOfGravity();
    // m_CentralMoments = src->GetCentralMoments();
    m_WeightedPrincipalMoments = src->GetWeightedPrincipalMoments();
    m_WeightedPrincipalAxes = src->GetWeightedPrincipalAxes();
    m_Kurtosis = src->GetKurtosis();
    m_Skewness = src->GetSkewness();
    m_WeightedElongation = src->GetWeightedElongation();
    m_Histogram = src->GetHistogram();
    m_WeightedFlatness = src->GetWeightedFlatness();
  }

  template <typename TSourceLabelObject>
  void
  CopyAllFrom(const TSourceLabelObject * src)
  {
    itkAssertOrThrowMacro((src != nullptr), "Null Pointer");
    this->template CopyLinesFrom<TSourceLabelObject>(src);
    this->template CopyAttributesFrom<TSourceLabelObject>(src);
  }

  const double &
  GetMinimum() const
  {
    return m_Minimum;
  }

  void
  SetMinimum(const double v)
  {
    m_Minimum = v;
  }

  const double &
  GetMaximum() const
  {
    return m_Maximum;
  }

  void
  SetMaximum(const double v)
  {
    m_Maximum = v;
  }

  const double &
  GetMean() const
  {
    return m_Mean;
  }

  void
  SetMean(const double v)
  {
    m_Mean = v;
  }

  const double &
  GetSum() const
  {
    return m_Sum;
  }

  void
  SetSum(const double v)
  {
    m_Sum = v;
  }

  const double &
  GetStandardDeviation() const
  {
    return m_StandardDeviation;
  }

  void
  SetStandardDeviation(const double v)
  {
    m_StandardDeviation = v;
  }

  const double &
  GetVariance() const
  {
    return m_Variance;
  }

  void
  SetVariance(const double v)
  {
    m_Variance = v;
  }

  const double &
  GetMedian() const
  {
    return m_Median;
  }

  void
  SetMedian(const double v)
  {
    m_Median = v;
  }

  const IndexType &
  GetMaximumIndex() const
  {
    return m_MaximumIndex;
  }

  void
  SetMaximumIndex(const IndexType & v)
  {
    m_MaximumIndex = v;
  }

  const IndexType &
  GetMinimumIndex() const
  {
    return m_MinimumIndex;
  }

  void
  SetMinimumIndex(const IndexType & v)
  {
    m_MinimumIndex = v;
  }

  const PointType &
  GetCenterOfGravity() const
  {
    return m_CenterOfGravity;
  }

  void
  SetCenterOfGravity(const PointType & v)
  {
    m_CenterOfGravity = v;
  }

  /*
  const MatrixType & GetCentralMoments() const
    {
    return m_CentralMoments;
    }

  void SetCentralMoments( const MatrixType & v )
    {
    m_CentralMoments = v;
    }*/

  const VectorType &
  GetWeightedPrincipalMoments() const
  {
    return m_WeightedPrincipalMoments;
  }

  void
  SetWeightedPrincipalMoments(const VectorType & v)
  {
    m_WeightedPrincipalMoments = v;
  }

  const MatrixType &
  GetWeightedPrincipalAxes() const
  {
    return m_WeightedPrincipalAxes;
  }

  void
  SetWeightedPrincipalAxes(const MatrixType & v)
  {
    m_WeightedPrincipalAxes = v;
  }

  const double &
  GetSkewness() const
  {
    return m_Skewness;
  }

  void
  SetSkewness(const double v)
  {
    m_Skewness = v;
  }

  const double &
  GetKurtosis() const
  {
    return m_Kurtosis;
  }

  void
  SetKurtosis(const double v)
  {
    m_Kurtosis = v;
  }

  const double &
  GetWeightedElongation() const
  {
    return m_WeightedElongation;
  }

  void
  SetWeightedElongation(const double v)
  {
    m_WeightedElongation = v;
  }

  const HistogramType *
  GetHistogram() const
  {
    return m_Histogram;
  }

  void
  SetHistogram(const HistogramType * v)
  {
    m_Histogram = v;
  }

  const double &
  GetWeightedFlatness() const
  {
    return m_WeightedFlatness;
  }

  void
  SetWeightedFlatness(const double v)
  {
    m_WeightedFlatness = v;
  }

  // some helper methods - not really required, but really useful!
  /** Affine transform for mapping to and from principal axis */
  using AffineTransformType = AffineTransform<double, Self::ImageDimension>;
  using AffineTransformPointer = typename AffineTransformType::Pointer;

  /** Get the affine transform from principal axes to physical axes
   * This method returns an affine transform which transforms from
   * the principal axes coordinate system to physical coordinates. */
  AffineTransformPointer
  GetWeightedPrincipalAxesToPhysicalAxesTransform() const
  {
    typename AffineTransformType::MatrixType matrix;
    typename AffineTransformType::OffsetType offset;
    for (unsigned int i = 0; i < ImageDimension; ++i)
    {
      offset[i] = m_CenterOfGravity[i];
      for (unsigned int j = 0; j < ImageDimension; ++j)
      {
        matrix[j][i] = m_WeightedPrincipalAxes[i][j]; // Note the transposition
      }
    }

    AffineTransformPointer result = AffineTransformType::New();

    result->SetMatrix(matrix);
    result->SetOffset(offset);

    return result;
  }

  /** Get the affine transform from physical axes to principal axes
   * This method returns an affine transform which transforms from
   * the physical coordinate system to the principal axes coordinate
   * system. */
  AffineTransformPointer
  GetPhysicalAxesToWeightedPrincipalAxesTransform() const
  {
    typename AffineTransformType::MatrixType matrix;
    typename AffineTransformType::OffsetType offset;
    for (unsigned int i = 0; i < ImageDimension; ++i)
    {
      offset[i] = m_CenterOfGravity[i];
      for (unsigned int j = 0; j < ImageDimension; ++j)
      {
        matrix[j][i] = m_WeightedPrincipalAxes[i][j]; // Note the transposition
      }
    }

    AffineTransformPointer result = AffineTransformType::New();
    result->SetMatrix(matrix);
    result->SetOffset(offset);

    AffineTransformPointer inverse = AffineTransformType::New();
    result->GetInverse(inverse);

    return inverse;
  }

protected:
  StatisticsLabelObject()
  {
    m_Minimum = 0;
    m_Maximum = 0;
    m_Mean = 0;
    m_Sum = 0;
    m_StandardDeviation = 0;
    m_Variance = 0;
    m_Median = 0;
    m_MaximumIndex.Fill(0);
    m_MinimumIndex.Fill(0);
    m_CenterOfGravity.Fill(0);
    // m_CentralMoments.Fill(0);
    m_WeightedPrincipalMoments.Fill(0);
    m_WeightedPrincipalAxes.Fill(0);
    m_Kurtosis = 0;
    m_Skewness = 0;
    m_WeightedElongation = 0;
    m_Histogram = nullptr;
    m_WeightedFlatness = 0;
  }

  void
  PrintSelf(std::ostream & os, Indent indent) const override
  {
    Superclass::PrintSelf(os, indent);

    os << indent << "Minimum: " << m_Minimum << std::endl;
    os << indent << "Maximum: " << m_Maximum << std::endl;
    os << indent << "Mean: " << m_Mean << std::endl;
    os << indent << "Sum: " << m_Sum << std::endl;
    os << indent << "StandardDeviation: " << m_StandardDeviation << std::endl;
    os << indent << "Variance: " << m_Variance << std::endl;
    os << indent << "Median: " << m_Median << std::endl;
    os << indent << "Skewness: " << m_Skewness << std::endl;
    os << indent << "Kurtosis: " << m_Kurtosis << std::endl;
    os << indent << "WeightedElongation: " << m_WeightedElongation << std::endl;
    os << indent << "WeightedFlatness: " << m_WeightedFlatness << std::endl;
    os << indent << "MaximumIndex: " << m_MaximumIndex << std::endl;
    os << indent << "MinimumIndex: " << m_MinimumIndex << std::endl;
    os << indent << "CenterOfGravity: " << m_CenterOfGravity << std::endl;
    // os << indent << "CentralMoments: " << std::endl << m_CentralMoments;
    os << indent << "WeightedPrincipalMoments: " << m_WeightedPrincipalMoments << std::endl;
    os << indent << "WeightedPrincipalAxes: " << std::endl << m_WeightedPrincipalAxes;
    itkPrintSelfObjectMacro(Histogram);
  }

private:
  double    m_Minimum{};
  double    m_Maximum{};
  double    m_Mean{};
  double    m_Sum{};
  double    m_StandardDeviation{};
  double    m_Variance{};
  double    m_Median{};
  IndexType m_MaximumIndex{};
  IndexType m_MinimumIndex{};
  PointType m_CenterOfGravity{};
  // MatrixType m_CentralMoments;
  VectorType m_WeightedPrincipalMoments{};
  MatrixType m_WeightedPrincipalAxes{};
  double     m_Skewness{};
  double     m_Kurtosis{};
  double     m_WeightedElongation{};

  typename HistogramType::ConstPointer m_Histogram{};

  double m_WeightedFlatness{};
};
} // end namespace itk

#endif