Commit fd1e8306 authored by Leon Merten Lohse's avatar Leon Merten Lohse

python version

parent dbe6a828
CXXFLAGS= -Wall -Wpedantic -std=c++14 -Iinclude/ -O3 -march=native -fopenmp
CXXFLAGS= -Wall -Wpedantic -std=c++14 -Iinclude/ -I/usr/include/eigen3/ -O3 -march=native -fopenmp
LDFLAGS=
LDLIBS= -lm -lstdc++ -ltiff -lgomp
LDLIBSPY= -lm -lstdc++ -lgomp
CXX=g++
EXS= moenchInterpolation moenchMakeEta histogram
......@@ -9,6 +10,11 @@ all: $(EXS)
HEADERS=$(wildcard include/*.h)
SRCS=moench03Interpolation.cpp
PY_INCLUDES=$(shell python3 -m pybind11 --includes)
PY_EXTENSION_SUFFIX=$(shell python3-config --extension-suffix)
PY_MODULE_NAME=_moench
moenchInterpolation: $(SRCS) $(HEADERS)
$(CXX) -o moenchInterpolation $(CXXFLAGS) $(SRCS) $(LDFLAGS) $(LDLIBS)
......@@ -18,5 +24,9 @@ moenchMakeEta: $(SRCS) $(HEADERS)
histogram: histogram.cpp $(HEADERS)
$(CXX) -o histogram $(CXXFLAGS) histogram.cpp $(LDFLAGS) $(LDLIBS)
python: python.cpp $(HEADERS)
$(CXX) $(CXXFLAGS) -shared -fPIC $(PY_INCLUDES) python.cpp -o $(PY_MODULE_NAME)$(PY_EXTENSION_SUFFIX) $(LDFLAGS) $(LDLIBSPY)
clean:
rm $(OBJS) $(EXS)
rm $(OBJS) $(EXS) $(PY_MODULE_NAME)$(PY_EXTENSION_SUFFIX)
......@@ -4,13 +4,16 @@
#include <vector>
#include <cassert>
#include <eigen3/Eigen/Dense>
#include <Eigen/Dense>
using CoordinateScalar = float;
using DataScalar = float;
using Coordinate = Eigen::Array<CoordinateScalar, 2, 1>;
#include "tiffIO.h"
template <typename T>
//using ImageData = Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
using ImageData = Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic>;
using ImageData = Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
template <typename T>
struct BaseImage {
......@@ -38,33 +41,6 @@ struct BaseImage {
return pixels(iy, ix);;
}
void write(const std::string & imgname) {
/* eigen stores 2d arrays column-major and seems to be very inefficient with row-major arrays
* TODO: find out why */
ImageData<T> tmp = pixels.transpose();
std::vector<float> gm (tmp.data(), tmp.data()+ tmp.size());
WriteToTiff(gm, pixels.cols(), pixels.rows(), imgname);
}
void read(const std::string & filename) {
auto tiffimage = ReadFromTiff(filename);
if (tiffimage) {
if (tiffimage->width == pixels.cols() && tiffimage->height == pixels.rows()) {
auto gm = tiffimage->data;
ImageData<T> tmp (tiffimage->width, tiffimage->height);
std::copy(std::begin(gm), std::end(gm), tmp.data());
pixels = tmp.transpose();
} else {
std::cout << "ERROR: image " << filename << " has wrong dimension: "
<< tiffimage->width << "x" << tiffimage->height
<< " instead of " << pixels.cols() << "x" << pixels.rows() << std::endl;
}
}else{
std::cout << "ERROR: image unreadable" << std::endl;
}
}
private:
ImageData<T> pixels;
};
......@@ -76,12 +52,12 @@ struct InterpolatedImage : BaseImage<T> {
InterpolatedImage(int nx, int ny, int ns)
: BaseImage<T>(nx * ns, ny * ns),
n_subpixels(ns)
n_subpixels(ns)
{ };
T& operator[](std::array<double, 2> position) {
int ix= static_cast<double>(n_subpixels) * position[0];
int iy= static_cast<double>(n_subpixels) * position[1];
T& operator[](std::array<CoordinateScalar, 2> position) {
int ix= static_cast<CoordinateScalar>(n_subpixels) * position[0];
int iy= static_cast<CoordinateScalar>(n_subpixels) * position[1];
return (*this)(ix, iy);
}
......@@ -89,18 +65,18 @@ struct InterpolatedImage : BaseImage<T> {
int n_subpixels;
};
template <int nb=400, int emin=0, int emax=1, typename T=int>
template <int emin=0, int emax=1, typename T=int>
struct Field : BaseImage<T> {
static constexpr int nbeta = nb;
static constexpr double etamin = emin;
static constexpr double etamax = emax;
static constexpr double etastep = (etamax-etamin)/nbeta;
static constexpr CoordinateScalar etamin = emin;
static constexpr CoordinateScalar etamax = emax;
static constexpr CoordinateScalar etastep = (etamax-etamin)/nbeta;
Field()
Field(int nb)
: BaseImage<T>(nb, nb) {}
inline static std::array<int, 2> etaToIndex(double etax, double etay){
inline static std::array<int, 2> etaToIndex(CoordinateScalar etax, CoordinateScalar etay){
assert(etax >= etamin);
assert(etax < etamax);
assert(etay >= etamin);
......@@ -110,25 +86,24 @@ struct Field : BaseImage<T> {
return { ix, iy };
}
T & operator[](std::array<double, 2> eta) {
T & operator[](std::array<CoordinateScalar, 2> eta) {
auto index = etaToIndex(eta[0], eta[1]);
return (*this)(index[0], index[1]); // (x,y)
}
};
template <int resolution=100>
struct Chart2d {
Field<resolution, 0, 1, double> e_x;
Field<resolution, 0, 1, double> e_y;
Field<0, 1, CoordinateScalar> e_x;
Field<0, 1, CoordinateScalar> e_y;
Chart2d(Field<resolution, 0, 1, double> _e_x, Field<resolution, 0, 1, double> _e_y)
Chart2d(Field<0, 1, CoordinateScalar> _e_x, Field<0, 1, CoordinateScalar> _e_y)
{
e_x = _e_x;
e_y = _e_y;
}
std::array<double,2> operator()(std::array<double, 2> chi) {
std::array<double,2> x;
std::array<CoordinateScalar,2> operator()(std::array<CoordinateScalar, 2> chi) {
std::array<CoordinateScalar,2> x;
x[0] = e_x[chi];
x[1] = e_y[chi];
......
This diff is collapsed.
......@@ -4,9 +4,25 @@
#include <array>
#include <iterator>
#include <algorithm>
#include <eigen3/Eigen/Dense>
#include <Eigen/Dense>
#include <iostream>
using DataScalar = float;
using ArrayXX = Eigen::Array<DataScalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
using ArrayX = Eigen::Array<DataScalar, Eigen::Dynamic, 1>;
using Array33 = Eigen::Array<DataScalar, 3, 3>;
using Array22 = Eigen::Array<DataScalar, 2, 2>;
using RawArray33 = Eigen::Array<int, 3, 3, Eigen::RowMajor>;
using gainmap = Eigen::Array<DataScalar, Eigen::Dynamic, Eigen::Dynamic>;
using CoordinateScalar = double;
using Coordinate = Eigen::Array<CoordinateScalar, 2, 1>;
template <int dim>
using Index = Eigen::Array<int, dim, 1>;
template <int dx=3, int dy=3>
struct single_photon_hit {
......@@ -40,26 +56,26 @@ int argmax(T container) {
/* Helper data structure to analyze single photon events */
struct photon_event {
double sum;
double totquad;
DataScalar sum;
DataScalar totquad;
std::array<double, 2> eta;
std::array<int, 2> position;
Coordinate eta;
Index<2> position;
photon_event(const single_photon_hit<3,3> &hit, const Eigen::ArrayXXd & gainmap) {
const Eigen::Map<Eigen::Array<int, 3, 3, Eigen::RowMajor>> raw_data (const_cast<int *>(&(hit.data[0])), 3, 3);
photon_event(const single_photon_hit<3,3> &hit, const Eigen::Ref<const gainmap> & gainmap) {
const Eigen::Map<RawArray33> raw_data (const_cast<int *>(&(hit.data[0])), 3, 3);
position = { hit.x - 1, hit.y - 1 };
Eigen::Array33d data = raw_data.cast<double>() * gainmap.block<3,3>(position[1], position[0]);
Array33 data = raw_data.cast<DataScalar>() * gainmap.block<3,3>(position[1], position[0]);
sum = data.sum();
// TODO: static machen
// offset of the 2x2-quadrant relative to the center pixel of the 3x3-cluster (x, y)
const std::array<std::array<int,2>,4> quadrant_offsets = {{ { 0, 0 }, { 0, 1}, { 1, 0 }, { 1, 1} }};
sum = data.sum();
const std::array<Index<2>,4> quadrant_offsets = {{ Index<2>(0, 0), Index<2>(0, 1), Index<2>(1, 0), Index<2>(1, 1) }};
std::array<double, 4> quadrant_sums;
std::array<DataScalar, 4> quadrant_sums;
// index in matrix convention (i, j) = (y,x)
quadrant_sums[0] = data.block<2,2>(0, 0).sum(); // top left
quadrant_sums[1] = data.block<2,2>(1, 0).sum(); // bottom left
......@@ -70,16 +86,15 @@ struct photon_event {
int maximal_quadrant = argmax(quadrant_sums);
totquad = quadrant_sums[maximal_quadrant];
std::array<int, 2> offset = quadrant_offsets[maximal_quadrant];
Index<2> offset = quadrant_offsets[maximal_quadrant];
const Eigen::Array22d sDum = data.block<2,2>(offset[1], offset[0]);
const Array22 sDum = data.block<2,2>(offset[1], offset[0]);
position[0] += offset[0]; // x
position[1] += offset[1]; // y
position += offset;
// calculate eta
eta[0]=(sDum(0,1) + sDum(1,1))/totquad; // x
eta[1]=(sDum(1,0) + sDum(1,1))/totquad; // y
eta(0)=(sDum(0,1) + sDum(1,1))/totquad; // x
eta(1)=(sDum(1,0) + sDum(1,1))/totquad; // y
}
};
......
#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>
#include <pybind11/stl.h>
#include "moench.h"
namespace py = pybind11;
PYBIND11_MODULE(_moench, m) {
m.doc() = "MOENCH detector io";
m.def("compute_hitmaps", &moench::compute_hitmaps, py::return_value_policy::move);
m.def("compute_image", &moench::compute_image, py::return_value_policy::move);
m.attr("dx") = py::int_(moench::dx);
m.attr("dy") = py::int_(moench::dy);
//m.attr("hitmap_resolution") = py::int_(moench::hitmap_resolution);
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment