2007-03-08 04:03:40 +08:00
|
|
|
/**********************************************************************
|
|
|
|
* File: linlsq.h (Formerly llsq.h)
|
|
|
|
* Description: Linear Least squares fitting code.
|
|
|
|
* Author: Ray Smith
|
|
|
|
* Created: Thu Sep 12 08:44:51 BST 1991
|
|
|
|
*
|
|
|
|
* (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.
|
|
|
|
*
|
|
|
|
**********************************************************************/
|
|
|
|
|
2011-03-22 05:44:45 +08:00
|
|
|
#ifndef TESSERACT_CCSTRUCT_LINLSQ_H_
|
|
|
|
#define TESSERACT_CCSTRUCT_LINLSQ_H_
|
2007-03-08 04:03:40 +08:00
|
|
|
|
2011-03-22 05:44:45 +08:00
|
|
|
#include "points.h"
|
|
|
|
#include "params.h"
|
2007-03-08 04:03:40 +08:00
|
|
|
|
2011-03-22 05:44:45 +08:00
|
|
|
class LLSQ {
|
|
|
|
public:
|
|
|
|
LLSQ() { // constructor
|
|
|
|
clear(); // set to zeros
|
|
|
|
}
|
|
|
|
void clear(); // initialize
|
2007-03-08 04:03:40 +08:00
|
|
|
|
2011-03-22 05:44:45 +08:00
|
|
|
// 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 count() const { // no of elements
|
|
|
|
return static_cast<int>(total_weight + 0.5);
|
|
|
|
}
|
2007-03-08 04:03:40 +08:00
|
|
|
|
2011-03-22 05:44:45 +08:00
|
|
|
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.
|
2007-03-08 04:03:40 +08:00
|
|
|
|
2011-03-22 05:44:45 +08:00
|
|
|
// Returns the x,y means as an FCOORD.
|
|
|
|
FCOORD mean_point() const;
|
2013-09-21 03:43:47 +08:00
|
|
|
|
|
|
|
// 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;
|
|
|
|
|
2011-03-22 05:44:45 +08:00
|
|
|
// 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
|
2013-09-21 03:43:47 +08:00
|
|
|
// 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.
|
2011-03-22 05:44:45 +08:00
|
|
|
FCOORD vector_fit() const;
|
2013-09-21 03:43:47 +08:00
|
|
|
|
2011-03-22 05:44:45 +08:00
|
|
|
// 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;
|
|
|
|
}
|
2007-03-08 04:03:40 +08:00
|
|
|
|
2011-03-22 05:44:45 +08:00
|
|
|
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
|
2007-03-08 04:03:40 +08:00
|
|
|
};
|
|
|
|
|
2013-09-21 03:43:47 +08:00
|
|
|
|
|
|
|
// 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 GenericVector, as it makes heavy used of LLSQ.
|
|
|
|
// T must be an integer or float/double type.
|
|
|
|
template<typename T> T MedianOfCircularValues(T modulus, GenericVector<T>* v) {
|
|
|
|
LLSQ stats;
|
|
|
|
T halfrange = static_cast<T>(modulus / 2);
|
|
|
|
int num_elements = v->size();
|
|
|
|
for (int i = 0; i < num_elements; ++i) {
|
|
|
|
stats.add((*v)[i], (*v)[i] + halfrange);
|
|
|
|
}
|
|
|
|
bool offset_needed = stats.y_variance() < stats.x_variance();
|
|
|
|
if (offset_needed) {
|
|
|
|
for (int i = 0; i < num_elements; ++i) {
|
|
|
|
(*v)[i] += halfrange;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
int median_index = v->choose_nth_item(num_elements / 2);
|
|
|
|
if (offset_needed) {
|
|
|
|
for (int i = 0; i < num_elements; ++i) {
|
|
|
|
(*v)[i] -= halfrange;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
return (*v)[median_index];
|
|
|
|
}
|
|
|
|
|
|
|
|
|
2011-03-22 05:44:45 +08:00
|
|
|
#endif // TESSERACT_CCSTRUCT_LINLSQ_H_
|