linefit.cpp

barcode readers [ from Image]

//
// Fit a line to a vector of data. Uses method from section 15.2 of
//   "Numerical Recipes in C : The Art of Scientific Computing"
//   by William H. Press, Brian P. Flannery, Saul A. Teukolsky, William T. Vetterling.
//   Cambridge University Press; 2 edition (October 30, 1992).
//   The original source of this code was the fit function from
//   Section 15.2 of the book. See http://library.lanl.gov/numerical/bookcpdf.html
//   for a complete PDF of the book.
//
// Copyright (C) 2006 by Jon A. Webb (Contact via GMail; username is jonawebb)
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// 
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
// 
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
//
//

#include "LineFit.h"

namespace Algorithm
{

    EXPORT_C CLineFit* CLineFit::NewL()
    {
        return new (ELeave) CLineFit();
    }

    CLineFit::CLineFit() 
    {
    }

    CLineFit::~CLineFit(void)
    {
    }

    bool CLineFit::FitExact(CArrayFix<TPoint>* pPoints, TReal *a, TReal *b, TReal *siga, TReal *sigb, TReal *chi2)
    {
		/*
			Given an array of TPoint's, 
			fit them to a straight line iY = a + b*iX by minimizing ?2. Returned are
			a,b and their respective probable uncertainties siga and sigb, the chi-square chi2, and the
			goodness-of-fit probability q (that the fit would have ?2 this large or larger). 
			The normalization of chi2 is to unit standard deviation on all points.
		*/
		int i;
		TReal t, sxoss, sx=0.0, sy=0.0, st2=0.0, ss, sigdat;
		*b=0.0;
		for (i=0; i<pPoints->Count(); i++) { // without weights.
			sx += pPoints->At(i).iX;
			sy += pPoints->At(i).iY;
		}
		ss = TReal(pPoints->Count());
		sxoss = sx/ss;
		for (i=0;i<pPoints->Count();i++) {
			t= pPoints->At(i).iX - sxoss;
			st2 += t*t;
			*b += t*pPoints->At(i).iY;
		}
		*b /= st2; // Solve for a, b, sa, and sb.
		*a=(sy-sx*(*b))/ss;
		Math::Sqrt(*siga, (1.0+sx*sx/(ss*st2))/ss);
		Math::Sqrt(*sigb, 1.0/st2);
		*chi2=0.0; // Calculate ?2.
		for (i=0; i<pPoints->Count(); i++) {
			TReal r = pPoints->At(i).iY-(*a)-(*b)*pPoints->At(i).iX;
			*chi2 += r*r;
		}
		Math::Sqrt(sigdat, (*chi2)/(pPoints->Count()-2)); // For unweighted data evaluate typical
		// sig using chi2, and adjust
		// the standard deviations.
		*siga *= sigdat;
		*sigb *= sigdat;
		return true;
    }

	bool CLineFit::FitNoisy(CArrayFix<TPoint>* pPoints, TReal *a, TReal *b, TReal maxChi2)
	{
		TReal siga, sigb, chi2;
        FitExact(pPoints, a, b, &siga, &sigb, &chi2);
		while (chi2 > maxChi2) {
			// find the point that differs most from the fit line
			int i;
			TReal maxErr = 0.0;
			int iWorst = -1;
			for (i=0; i<pPoints->Count(); i++) {
				TReal err = *a * pPoints->At(i).iX + *b - pPoints->At(i).iY;
				if (err < 0) {
					err = -err;
				}
				if (err > maxErr) {
					iWorst = i;
					maxErr = err;
				}
			}
			// remove this point from the array
			pPoints->Delete(iWorst);
			FitExact(pPoints, a, b, &siga, &sigb, &chi2);
		}
		return true;
	}

};