Merge pull request #3 from hz-b/dev/feature/volume-add-save-cast

[TASK] volume of different integer types with internal range check

examples/plot_simple.py

@@ -2,39 +2,40 @@
 from pathlib import Path
 import numpy as np
 import matplotlib.pyplot as plt
-from timepix_sort.config import TDC2TriggerMode
-from timepix_sort.data_model import PixelEvent
-from timepix_sort.read import read_chunks
-from timepix_sort.process_chunks import process_chunks
+from timepix_sort.config import TDC2TriggerMode, TDC1TriggerMode
+from timepix_sort.read import read_chunks, read_file_to_buffer
+from timepix_sort.process import process_chunks
+from timepix_sort.post_process import data_to_volume
 
 data_dir = Path(__file__).parent.parent / "tests" / "data"
 filename = data_dir / "Co_pos_0000.tpx3"
 # filename = data_dir / "Co_test_0000.tpx3"
 
 
 start = datetime.datetime.now()
-events = [
-    ev
-    for ev in process_chunks(
-        read_chunks(filename), trigger_mode=TDC2TriggerMode.rising_edge, tot_min=1
-    )
-    if ev is not None
-]
+chunks = read_chunks(read_file_to_buffer(filename))
+events, event_statistics = process_chunks(chunks, TDC1TriggerMode.rising_edge, 6)
 end = datetime.datetime.now()
 dt = (end - start).total_seconds()
 print(f"Processing {len(events)} events required {dt:.3f} seconds")
 
-result = np.zeros([515, 514])
-for ev in events:
-    if isinstance(ev, PixelEvent):
-        x, y = ev.pos.x, ev.pos.y
-        assert x >= 0
-        assert y >= 0
-        x = int(round(x))
-        y = int(round(y))
-        result[x, y] += 1
+si_buf = events.sorted_indices()
+sorted_indices = np.array(si_buf, copy=False)
 
+pixels_diff = events.pixel_events_with_difference_time(si_buf)
+pixels_diff.sort()
+assert pixels_diff.is_sorted
+
+tmp = np.arange(0, 20)
+lut = np.array([tmp, tmp]).T
+lut[:, 0] *= int(1.6e6 / 20)
+volume = np.zeros([525, 524, len(tmp)], dtype=np.int16)
+data_to_volume(pixels_diff, lut, volume)
+result = np.sum(volume, axis=-1)
+largest_hit = np.max(result.ravel())
+lowest_hit = np.min(result.ravel())
 np.save("image_data.npy", result)
-plt.imshow(result[100:350, :400])
+print(f"hit range {lowest_hit} {largest_hit}")
+plt.imshow(result[100:350, :400], vmax=largest_hit/3)
 plt.axis("equal")
 plt.savefig("rough_data_result.png")

examples/use_ext.py

@@ -44,9 +44,9 @@
     events_statistics.n_control_indications
 )
 txt = f"""Event statistics
-# of events  {events_statistics.n_events: 8d} 
+# of events  {events_statistics.n_events: 8d}
 
-consisting of 
+consisting of
 pixel        {events_statistics.n_pixels: 8d}
 timestamps   {events_statistics.n_timestamps: 8d}
 global time  {events_statistics.n_global_time: 8d}
@@ -70,9 +70,9 @@
 toc_pixels_sorted = _now()
 
 tmp = np.arange(0, 20)
-lut = np.array([tmp, tmp]).T
+lut = np.array([tmp, tmp], dtype=np.uint64).T
 lut[:, 0] *= int(1.6e6 / 20)
-volume = np.zeros([525, 524, len(tmp)], dtype=np.uint16)
+volume = np.zeros([525, 524, len(tmp)], dtype=np.int32)
 data_to_volume(pixels_diff, lut, volume)
 # print("volume sum", np.sum(np.sum(volume)))
 

src/c++/volume.cpp

@@ -3,16 +3,19 @@
 #include <pybind11/numpy.h>
 #include <timepix_sort/events.h>
 #include <iostream>
+#include <boost/numeric/conversion/cast.hpp>
 
 namespace py = pybind11;
 namespace dm = timepix::data_model;
 namespace tpp = timepix::python;
 namespace tps = timepix::sort;
 
 
+template<typename array_base_t>
 static void
-data_to_volume(const tps::PixelEventsDiffTime& pd, const py::array_t<uint64_t>& lut, py::array_t<uint16_t> volume)
+data_to_volume(const tps::PixelEventsDiffTime& pd, const py::array_t<uint64_t>& lut, py::array& volume)
 {
+
     auto lut_r = lut.unchecked<2>();
     if (lut_r.shape(1) != 2)
 	throw std::runtime_error("Lut last dimension must be 2");
@@ -24,9 +27,11 @@ data_to_volume(const tps::PixelEventsDiffTime& pd, const py::array_t<uint64_t>&
 	lutt.insert(lut_r(i, 0), lut_r(i, 1));
     }
 
-    auto r = volume.mutable_unchecked<3>();
+    auto r = volume.mutable_unchecked<array_base_t, 3>();
     // fill the buffer
+    size_t event_nr = -1;
     for(const auto& ev: pd){
+	event_nr++;
 	const auto pos = tps::map_pixel_and_chip_nr_to_global_pixel(ev.pos(), ev.chip_nr());
 	if (pos.x() < 0){
 	    throw std::range_error("x < 0");
@@ -40,7 +45,7 @@ data_to_volume(const tps::PixelEventsDiffTime& pd, const py::array_t<uint64_t>&
 	}
 	if (pos.y() >= r.shape(1) - 1){
 	    throw std::range_error("y >= r.shape(1) -1");
-		}
+	}
 	const auto tmp = lutt.linear_interp(ev.time_of_arrival());
 	const auto ti = int64_t(std::round(tmp));
 	if (ti < 0){
@@ -49,14 +54,36 @@ data_to_volume(const tps::PixelEventsDiffTime& pd, const py::array_t<uint64_t>&
 	if (ti >= r.shape(2) - 1){
 	    throw std::range_error("ti >= r.shape(2) - 1");
 	}
-	r(pos.x(), pos.y(), ti) += ev.time_over_threshold();
+	const int64_t tot = ev.time_over_threshold();
+	const int64_t prev_tot = r(pos.x(), pos.y(), ti);
+	const int64_t tot_s = tot + prev_tot;
+
+	try {
+	    r(pos.x(), pos.y(), ti) = boost::numeric_cast<array_base_t>(tot_s);
+	}
+	catch(boost::numeric::negative_overflow& e) {
+	    std::stringstream strm;
+	    strm << "event nr " << event_nr << ": tot = " << tot
+		 << " old value =  " << prev_tot
+		 << " total " << tot_s
+		 << ": " << e.what();
+	    throw std::overflow_error(strm.str());
+	}
+	catch(boost::numeric::positive_overflow& e) {
+	    std::stringstream strm;
+	    strm << "event nr " << event_nr << ": tot = " << tot
+		 << " old value =  " << prev_tot
+		 << " total " << tot_s
+		 << ": " << e.what();
+	    throw std::overflow_error(strm.str());
+	}
     }
 }
 
 static auto
 data_to_points(const tps::PixelEventsDiffTime& pd)
 {
-    std::vector<py::ssize_t> dims = {pd.size(), 3};
+    std::vector<py::ssize_t> dims = {boost::numeric_cast<py::ssize_t>(pd.size()), 3};
     py::array_t<uint64_t> res(dims);
     auto r = res.mutable_unchecked<2>();
 
@@ -75,9 +102,41 @@ data_to_points(const tps::PixelEventsDiffTime& pd)
 }
 
 
+static void
+data_to_volume_gateway(const tps::PixelEventsDiffTime& pd, const py::array_t<uint64_t>& lut, py::array& volume)
+{
+
+    const auto dtype = volume.dtype();
+    const auto dtype_num = dtype.num();
+    if(dtype_num == py::dtype::of<uint8_t>().num()){
+	data_to_volume<uint8_t>(pd, lut, volume);
+    } else if (dtype_num == py::dtype::of<int8_t>().num()){
+	data_to_volume<int8_t>(pd, lut, volume);
+    } else if(dtype_num == py::dtype::of<uint16_t>().num()){
+	data_to_volume<uint16_t>(pd, lut, volume);
+    } else if (dtype_num == py::dtype::of<int16_t>().num()){
+	data_to_volume<int16_t>(pd, lut, volume);
+    } else if(dtype_num == py::dtype::of<uint32_t>().num()){
+	data_to_volume<uint32_t>(pd, lut, volume);
+    } else if (dtype_num == py::dtype::of<int32_t>().num()){
+	data_to_volume<int32_t>(pd, lut, volume);
+    } else if(dtype_num ==  py::dtype::of<uint64_t>().num()){
+	data_to_volume<uint64_t>(pd, lut, volume);
+    } else if (dtype_num == py::dtype::of<int64_t>().num()){
+	data_to_volume<int64_t>(pd, lut, volume);
+    } else {
+	std::stringstream strm;
+	strm <<  "data to gateway dtype: " << dtype.char_()
+	     << ", dtype num " << dtype.num()
+	     << "not handled!";
+	throw std::runtime_error(strm.str());
+    }
+}
+
+
 void tpp::volume_init(pybind11::module &m)
 {
 
-    m.def("data_to_volume", data_to_volume);
+    m.def("data_to_volume", &data_to_volume_gateway);
     m.def("data_to_points", data_to_points);
 }