/**********************************************************************

 * File:        linlsq.h  (Formerly llsq.h)

 * Description: Linear Least squares fitting code.

 * Author:      Ray Smith

 *

 * (C) Copyright 1991, Hewlett-Packard Ltd.

 ** 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

 ** http://www.apache.org/licenses/LICENSE-2.0

 ** 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 TESSERACT_CCSTRUCT_LINLSQ_H_

#define TESSERACT_CCSTRUCT_LINLSQ_H_

#include "points.h" // for FCOORD

#include <algorithm> // for std::nth_element

#include <cstdint> // for int32_t

namespace tesseract {

class TESS_API LLSQ {

public:

  LLSQ() {   // constructor

    clear(); // set to zeros

  }

  void clear(); // initialize

  // Adds an element with a weight of 1.

  void add(double x, double y);

  // Adds an element with a specified weight.

  void add(double x, double y, double weight);

  // Adds a whole LLSQ.

  void add(const LLSQ &other);

  // Deletes an element with a weight of 1.

  void remove(double x, double y);

  int32_t count() const { // no of elements

    return static_cast<int>(total_weight + 0.5);

  }

  double m() const;                     // get gradient

  double c(double m) const;             // get constant

  double rms(double m, double c) const; // get error

  double pearson() const;               // get correlation coefficient.

  // Returns the x,y means as an FCOORD.

  FCOORD mean_point() const;

  // Returns the average sum of squared perpendicular error from a line

  // through mean_point() in the direction dir.

  double rms_orth(const FCOORD &dir) const;

  // Returns the direction of the fitted line as a unit vector, using the

  // least mean squared perpendicular distance. The line runs through the

  // mean_point, i.e. a point p on the line is given by:

  // p = mean_point() + lambda * vector_fit() for some real number lambda.

  // Note that the result (0<=x<=1, -1<=y<=-1) is directionally ambiguous

  // and may be negated without changing its meaning, since a line is only

  // unique to a range of pi radians.

  // Modernists prefer to think of this as an Eigenvalue problem, but

  // Pearson had the simple solution in 1901.

  //

  // Note that this is equivalent to returning the Principal Component in PCA,

  // or the eigenvector corresponding to the largest eigenvalue in the

  // covariance matrix.

  FCOORD vector_fit() const;

  // Returns the covariance.

  double covariance() const {

    if (total_weight > 0.0) {

      return (sigxy - sigx * sigy / total_weight) / total_weight;

    } else {

      return 0.0;

    }

  }

  double x_variance() const {

    if (total_weight > 0.0) {

      return (sigxx - sigx * sigx / total_weight) / total_weight;

    } else {

      return 0.0;

    }

  }

  double y_variance() const {

    if (total_weight > 0.0) {

      return (sigyy - sigy * sigy / total_weight) / total_weight;

    } else {

      return 0.0;

    }

  }

private:

  double total_weight; // no of elements or sum of weights.

  double sigx;         // sum of x

  double sigy;         // sum of y

  double sigxx;        // sum x squared

  double sigxy;        // sum of xy

  double sigyy;        // sum y squared

};

// Returns the median value of the vector, given that the values are

// circular, with the given modulus. Values may be signed or unsigned,

// eg range from -pi to pi (modulus 2pi) or from 0 to 2pi (modulus 2pi).

// NOTE that the array is shuffled, but the time taken is linear.

// An assumption is made that most of the values are spread over no more than

// half the range, but wrap-around is accounted for if the median is near

// the wrap-around point.

// Cannot be a member of vector, as it makes heavy use of LLSQ.

// T must be an integer or float/double type.

template <typename T>

T MedianOfCircularValues(T modulus, std::vector<T> &v) {

  LLSQ stats;

  T halfrange = static_cast<T>(modulus / 2);

  auto num_elements = v.size();

  for (auto i : v) {

    stats.add(i, i + halfrange);

  }

  bool offset_needed = stats.y_variance() < stats.x_variance();

  if (offset_needed) {

    for (auto i : v) {

      i += halfrange;

    }

  }

  auto median_index = num_elements / 2;

  std::nth_element(v.begin(), v.begin() + median_index, v.end());

  if (offset_needed) {

    for (auto i : v) {

      i -= halfrange;

    }

  }

  return v[median_index];

}

} // namespace tesseract

#endif // TESSERACT_CCSTRUCT_LINLSQ_H_