Add XYStat for Linear Regression calculation

include/SlsDetectorDefs.h

@@ -408,8 +408,35 @@ private:
 	static bool DoHist;
 };
  
+struct XYStat {
+	double xacc, xacc2, yacc, xyacc;
+	int xn;
+	double factor;
+	mutable Mutex lock;
+
+	XYStat(double f = 1);
+	void reset();
+	void add(double x, double y);
+	XYStat& operator =(const XYStat& o);
+	XYStat& operator += (const XYStat& o);
+
+	// Linear Regression
+	struct LinRegress {
+		int n;
+		double slope;
+		double offset;
+
+		LinRegress() : n(0), slope(0), offset(0)
+		{}
+	};
+
+	int n() const;
+	LinRegress calcLinRegress() const;
+};
+
 std::ostream& operator <<(std::ostream& os, const SimpleStat::Histogram& s);
 std::ostream& operator <<(std::ostream& os, const SimpleStat& s);
+std::ostream& operator <<(std::ostream& os, const XYStat::LinRegress& r);
 
 
 class Camera;

src/SlsDetectorDefs.cpp

@@ -475,6 +475,76 @@ double SimpleStat::std() const
 	return xn ? sqrt(xacc2 / xn - pow(ave(), 2)) : 0; 
 }
 
+XYStat::XYStat(double f)
+	: factor(f)
+{
+	reset();
+}
+
+void XYStat::reset()
+{
+	AutoMutex l(lock);
+	xacc = xacc2 = yacc = xyacc = 0;
+	xn = 0;
+}
+
+void XYStat::add(double x, double y) {
+	AutoMutex l(lock);
+	y *= factor;
+	xacc += x;
+	xacc2 += pow(x, 2);
+	yacc += y;
+	xyacc += x * y;
+	++xn;
+}
+
+XYStat& XYStat::operator =(const XYStat& o)
+{
+	if (&o == this)
+		return *this;
+
+	AutoMutex l(o.lock);
+	xacc = o.xacc;
+	xacc2 = o.xacc2;
+	yacc = o.yacc;
+	xyacc = o.xyacc;
+	xn = o.xn;
+	factor = o.factor;
+	return *this;
+}
+
+XYStat& XYStat::operator +=(const XYStat& o)
+{
+	if (o.factor != factor)
+		throw LIMA_HW_EXC(Error, "Cannot add different XYStats");
+
+	AutoMutex l(o.lock);
+	xacc += o.xacc;
+	xacc2 += o.xacc2;
+	yacc += o.yacc;
+	xyacc += o.xyacc;
+	xn += o.xn;
+	return *this;
+}
+
+int XYStat::n() const
+{ 
+	AutoMutex l(lock);
+	return xn;
+}
+
+XYStat::LinRegress XYStat::calcLinRegress() const
+{ 
+	AutoMutex l(lock);
+	LinRegress r;
+	if (!xn)
+		return r;
+	r.n = xn;
+	r.slope = (xn * xyacc - xacc * yacc) / (xn * xacc2 - xacc * xacc);
+	r.offset = (yacc - r.slope * xacc) / xn;
+	return r;
+}
+
 FrameType lima::SlsDetector::getLatestFrame(const FrameArray& l)
 {
 	DEB_STATIC_FUNCT();
@@ -569,6 +639,16 @@ ostream& lima::SlsDetector::operator <<(ostream& os, const SimpleStat& s)
 	return os << ">";
 }
 
+ostream& lima::SlsDetector::operator <<(ostream& os,
+					const XYStat::LinRegress& r)
+{
+	os << "<";
+	os << "slope=" << r.slope << ", "
+	   << "offset=" << r.offset << ", "
+	   << "n=" << r.n;
+	return os << ">";
+}
+
 TimeRangesChangedCallback::TimeRangesChangedCallback()
 	: m_cam(NULL)
 {

test/CMakeLists.txt

@@ -23,7 +23,8 @@
 set(test_src test_slsdetector 
              test_slsdetector_control
              test_thread_cpu_affinity
-	     test_expand_4)
+	     test_expand_4
+	     test_stat)
 
 limatools_run_camera_tests("${test_src}" ${PROJECT_NAME})
 

test/test_stat.cpp

+//###########################################################################
+// This file is part of LImA, a Library for Image Acquisition
+//
+// Copyright (C) : 2009-2011
+// European Synchrotron Radiation Facility
+// BP 220, Grenoble 38043
+// FRANCE
+//
+// This is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 3 of the License, or
+// (at your option) any later version.
+//
+// This software 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 General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, see <http://www.gnu.org/licenses/>.
+//###########################################################################
+
+#include "SlsDetectorDefs.h"
+#include "lima/Exceptions.h"
+#include <cmath>
+
+DEB_GLOBAL(DebModTest);
+
+using namespace std;
+using namespace lima;
+using namespace SlsDetector;
+
+bool is_approx(double x, double y, double tol)
+{
+	return fabs(x - y) <= tol;
+}
+
+bool is_approx_rel(double x, double y, double rel_tol)
+{
+	return is_approx(x, y, fabs(rel_tol * (fabs(x) + fabs(y)) / 2));
+}
+
+void test_stat(double slope, double offset, int n)
+{
+	DEB_GLOBAL_FUNCT();
+
+	XYStat stat;
+	for (int i = 0; i < n; ++i)
+		stat.add(i, slope * i + offset);
+	XYStat::LinRegress r = stat.calcLinRegress();
+	DEB_ALWAYS() << DEB_VAR1(r);
+	if (r.n != n)
+		THROW_HW_ERROR(Error) << DEB_VAR2(r.n, n);
+	else if (!is_approx_rel(r.slope, slope, 1e-3))
+		THROW_HW_ERROR(Error) << DEB_VAR2(r.slope, slope);
+	else if (!is_approx_rel(r.offset, offset, 1e-3))
+		THROW_HW_ERROR(Error) << DEB_VAR2(r.offset, offset);
+}
+
+int main(int argc, char *argv[])
+{
+	DEB_GLOBAL_FUNCT();
+
+	double slope, offset;
+	int n;
+
+	slope = 0.12345;
+	offset = 123.45;
+	n = 100;
+	test_stat(slope, offset, n);
+
+	return 0;
+}