From 506d3479498da9c56370915b421dcdffc85f2897 Mon Sep 17 00:00:00 2001
From: PWhiddy <all.cows.like.to.moo@gmail.com>
Date: Tue, 12 Sep 2017 16:46:38 -0700
Subject: [PATCH] PSF array can be set with set_array()

---
 search/src/PointSpreadFunc.cpp | 79 ++++++++++++++++++++++++++--------
 search/src/PointSpreadFunc.h   | 17 +++++++-
 2 files changed, 77 insertions(+), 19 deletions(-)

diff --git a/search/src/PointSpreadFunc.cpp b/search/src/PointSpreadFunc.cpp
index c91d77ddb..c267c077d 100644
--- a/search/src/PointSpreadFunc.cpp
+++ b/search/src/PointSpreadFunc.cpp
@@ -31,46 +31,91 @@ PointSpreadFunc::PointSpreadFunc(float stdev) {
 		i++;
 	}
 
-	radius = i;
-	dim = 2 * radius - 1;
-	sum = 0.0;
+	radius = i-1; // This value is good for
+	dim = 2*radius+1;
 
 	// Create 2D gaussain by multiplying with itself
 	kernel = std::vector<float>();
 	for (int ii = 0; ii < dim; ++ii) {
 		for (int jj = 0; jj < dim; ++jj) {
-			float current = simpleGauss[abs(radius - ii - 1)]
-					* simpleGauss[abs(radius - jj - 1)];
+			float current = simpleGauss[abs(radius - ii)]
+					* simpleGauss[abs(radius - jj)];
 			kernel.push_back(current);
-			sum += current;
 		}
 	}
+	calcSum();
+}
 
+PointSpreadFunc::PointSpreadFunc(PointSpreadFunc& other)
+{
+	kernel = other.getKernel();
+	dim = other.getDim();
+	radius = other.getRadius();
+	width = other.getStdev();
+	calcSum();
 }
 
-void PointSpreadFunc::squarePSF() {
+#ifdef Py_PYTHON_H
+PointSpreadFunc::PointSpreadFunc(pybind11::array_t<float> arr)
+{
+	setArray(arr);
+}
+
+void PointSpreadFunc::setArray(pybind11::array_t<float> arr)
+{
+	pybind11::buffer_info info = arr.request();
+
+	if (info.ndim != 2)
+		throw std::runtime_error("Array must have 2 dimensions. (It "
+				" must also be a square with odd dimensions)");
+
+	if (info.shape[0] != info.shape[1])
+		throw std::runtime_error(
+			"Array must be square (x-dimension == y-dimension)."
+			"It also must have odd dimensions.");
+	float *pix = static_cast<float*>(info.ptr);
+	dim = info.shape[0];
+	if (dim%2==0)
+		throw std::runtime_error("Array dimension must be odd. The "
+				"middle of an even numbered array is ambiguous.");
+	radius = dim/2; // Rounds down
 	sum = 0.0;
+	kernel = std::vector<float>(pix,pix+dim*dim);
+	calcSum();
+	width = 0.0;
+}
+#endif
+
+void PointSpreadFunc::calcSum()
+{
+	sum = 0.0;
+	for (auto& i : kernel) sum += i;
+}
+
+void PointSpreadFunc::squarePSF() {
 	for (float& i : kernel)
 	{
 		i = i * i;
-		sum += i;
 	}
+	calcSum();
 }
 
-void PointSpreadFunc::printPSF() {
-	std::cout.setf(std::ios::fixed, std::ios::floatfield);
-	std::cout.precision(3);
+std::string PointSpreadFunc::printPSF() {
+	std::stringstream ss;
+	ss.setf(std::ios::fixed, std::ios::floatfield);
+	ss.precision(3);
 	for (int row = 0; row < dim; ++row) {
-		std::cout << "| ";
+		ss << "| ";
 		for (int col = 0; col < dim; ++col) {
-			std::cout << kernel[row * dim + col] << " | ";
+			ss << kernel[row * dim + col] << " | ";
 		}
-		std::cout << "\n ";
+		ss << "\n ";
 		for (int space = 0; space < dim * 8 - 1; ++space)
-			std::cout << "-";
-		std::cout << "\n";
+			ss << "-";
+		ss << "\n";
 	}
-	std::cout << 100.0 * sum << "% of PSF contained within kernel\n";
+	ss << 100.0 * sum << "% of PSF contained within kernel\n";
+	return ss.str();
 }
 
 } /* namespace kbmod */
diff --git a/search/src/PointSpreadFunc.h b/search/src/PointSpreadFunc.h
index 7aebd7ff0..545df8933 100644
--- a/search/src/PointSpreadFunc.h
+++ b/search/src/PointSpreadFunc.h
@@ -12,6 +12,12 @@
 #include <iostream>
 #include <sstream>
 #include <vector>
+#include <stdexcept>
+#ifdef Py_PYTHON_H
+#include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
+#include <pybind11/stl.h>
+#endif
 #include "common.h"
 
 namespace kbmod {
@@ -20,16 +26,23 @@ class PointSpreadFunc
 {
 	public:
 		PointSpreadFunc(float stdev);
+		PointSpreadFunc(PointSpreadFunc& other);
+#ifdef Py_PYTHON_H
+		PointSpreadFunc(pybind11::array_t<float> arr);
+		void setArray(pybind11::array_t<float> arr);
+#endif
 		virtual ~PointSpreadFunc() {};
 		float getStdev() { return width; }
+		void calcSum();
 		float getSum() { return sum; }
 		int getDim() { return dim; }
-		int getRadius() { return radius-1; }
+		int getRadius() { return radius; }
 		int getSize() { return kernel.size(); }
 		std::vector<float> getKernel() { return kernel; };
 		float* kernelData() { return kernel.data(); }
 		void squarePSF();
-		void printPSF();
+		std::string printPSF();
+		// void normalize(); ???
 	private:
 		std::vector<float> kernel;
 		float width;