Update to new geometry system

CMakeLists.txt

@@ -56,7 +56,7 @@ set(
     "read_tiff.cpp"
     "dose_deposition.cpp"
     "geometry_benchmark.cpp"
-    "partitioning_tests.cpp"
+    #"partitioning_tests.cpp"
 )
 
 set(
@@ -90,8 +90,8 @@ if(WITH_ZMQ)
     add_executable("read_geometry_test" "examples/read_geometry_test.cpp")
     target_link_libraries("read_geometry_test" ${LIB_NAMES})
 
-    add_executable("read_legoscan_test" "examples/read_legoscan_test.cpp")
-    target_link_libraries("read_legoscan_test" ${LIB_NAMES})
+    #    add_executable("read_legoscan_test" "examples/read_legoscan_test.cpp")
+    #    target_link_libraries("read_legoscan_test" ${LIB_NAMES})
 
     add_executable("partition" "examples/partition.cpp")
     target_link_libraries("partition" ${LIB_NAMES})

distributed/distributed_reconstruction.cpp

@@ -35,7 +35,7 @@ void run(tomo::util::args opt) {
 
         // projectors and geometries are modified so that they are
         // intersected with volumes at proper location
-        auto geom = tomo::geometry::parallel<D, T>(global_volume, opt.k, opt.k);
+        auto geom = tomo::geometry::parallel<D, T>(global_volume, opt.k);
 
         // set up the partitioning
         std::array<int, D> size{};

distributed/simple_distributed.cpp

@@ -21,10 +21,10 @@ auto pixels(shadow s) {
     return (s.max_pt.x - s.min_pt.x) * (s.max_pt.y - s.min_pt.y);
 }
 
-math::vec<2_D, int> compute_pixel_intersection(math::vec3<T> point,
+math::vec<2_D, int> compute_pixel_intersection(math::vec3<T> x,
                                                geometry::projection<3_D, T> p) {
     auto geometric_point = math::ray_plane_intersection<T>(
-        p.source_location, point - p.source_location, p.detector_location,
+        p.source_location, x - p.source_location, p.detector_location,
         math::cross<T>(p.detector_tilt[0], p.detector_tilt[1]));
     assert(geometric_point);
     auto pt = geometric_point.value();
@@ -38,64 +38,61 @@ math::vec<2_D, int> compute_pixel_intersection(math::vec3<T> point,
             (int)(dy / p.detector_size[1] * p.detector_shape[1])};
 }
 
-shadow compute_shadow(std::vector<math::vec3<T>> points,
-                      geometry::projection<3_D, T> p, math::vec2<int> shape) {
+shadow compute_shadow(std::vector<math::vec3<T>> xs,
+                      geometry::projection<3_D, T> p) {
     shadow s = {
         {std::numeric_limits<int>::max(), std::numeric_limits<int>::max()},
         {-1, -1}};
-    for (auto pt : points) {
-        auto intersection = compute_pixel_intersection(pt, p);
-        intersection = math::min(intersection, shape - math::vec2<int>{1, 1});
+    for (auto x : xs) {
+        auto intersection = compute_pixel_intersection(x, p);
+        intersection =
+            math::min(intersection, p.detector_shape - math::vec2<int>{1, 1});
         intersection = math::max(intersection, {0, 0});
         s.min_pt = math::min(s.min_pt, intersection);
         s.max_pt = math::max(s.max_pt, intersection);
     }
+    // FIXME enlarge by one
     return s;
 }
 
 class restricted_geometry : public geometry::base<3, T> {
   public:
     restricted_geometry(geometry::trajectory<3_D, T>& geometry,
                         volume<3_D, T> local_volume)
-        : geometry::base<3, T>(geometry.projection_count(), {-1, -1}),
+        : geometry::base<3, T>(geometry.projection_count(), false),
           geometry_(geometry), local_volume_(local_volume),
           pixels_(geometry.projection_count() + 1) {
         project();
         std::transform(shadows_.begin(), shadows_.end(), pixels_.begin() + 1,
                        [&](auto shadow) { return pixels(shadow); });
         std::partial_sum(pixels_.begin(), pixels_.end(), pixels_.begin());
-        this->set_line_count(pixels_[pixels_.size() - 1]);
+        this->compute_lines_();
     }
 
     void project() {
         for (int i = 0; i < geometry_.projection_count(); ++i) {
             shadows_.push_back(compute_shadow(local_volume_.corners(),
-                                              geometry_.get_projection(i),
-                                              geometry_.detector_shape()));
+                                              geometry_.get_projection(i)));
         }
     }
 
-    virtual math::ray<3, T> get_line(int i) const override {
-        // 1 compute projection (FIXME how to identify in O(1)?, need to rethink
-        // the whole 'get_line' system, making the iterators locally can boost
-        // performance greatly)
-        auto p = (std::find_if(pixels_.begin(), pixels_.end(),
-                               [=](auto x) { return x > i; }) -
-                  pixels_.begin()) -
-                 1;
-        assert(p >= 0);
-        assert(p < geometry_.projection_count());
-        // 2 compute global index in projection
-        auto j = i - pixels_[p];
-        auto x = j % (shadows_[p].max_pt.x - shadows_[p].min_pt.x);
-        auto y = j / (shadows_[p].max_pt.x - shadows_[p].min_pt.x);
-        x += shadows_[p].min_pt.x;
-        y += shadows_[p].min_pt.y;
-        int global_index = p * geometry_.detector_pixel_count() +
-                           y * geometry_.detector_shape()[0] + x;
-        assert(global_index >= 0 && global_index < geometry_.lines());
-        // 3 return get_line of global index
-        return geometry_.get_line(global_index);
+    math::vec<2_D, int> projection_shape(int i) const override {
+        return {shadows_[i].max_pt.x - shadows_[i].min_pt.x,
+                shadows_[i].max_pt.y - shadows_[i].min_pt.y};
+    }
+
+    math::vec<3_D, T> detector_corner(int i) const override {
+        return geometry_.detector_corner(i) +
+               (T)shadows_[i].min_pt[0] * geometry_.projection_delta(i)[0] +
+               (T)shadows_[i].min_pt[1] * geometry_.projection_delta(i)[1];
+    }
+
+    math::vec<3_D, T> source_location(int i) const override {
+        return geometry_.source_location(i);
+    }
+
+    std::array<math::vec<3_D, T>, 2> projection_delta(int i) const override {
+        return geometry_.projection_delta(i);
     }
 
   private:
@@ -106,7 +103,8 @@ class restricted_geometry : public geometry::base<3, T> {
 };
 
 // Assumes global volume has physical size [0, 1]^3
-auto calculate_local_volume(bulk::rectangular_partitioning<3, 1>& partitioning,
+template <int G>
+auto calculate_local_volume(bulk::rectangular_partitioning<3, G>& partitioning,
                             int s) {
     auto to_vec = [](auto array_like) -> math::vec3<T> {
         math::vec3<T> result;
@@ -136,8 +134,8 @@ void run(util::args options) {
                                {options.k, options.k}, 10.0, 2.0);
 
     auto processors = env.available_processors();
-    auto partitioning = bulk::block_partitioning<3, 1>(
-        {options.k, options.k, options.k}, {processors}, {1});
+    auto partitioning = bulk::block_partitioning<3, 3>(
+        {options.k, options.k, options.k}, {2, 2, 2});
 
     env.spawn(processors, [&](auto& world) {
         auto local_volume =
@@ -170,6 +168,8 @@ void run(util::args options) {
         auto sinos = bulk::gather_all(rsino.data()); */
 
         plotter.plot(phantom);
+        plotter.send_projection_data(local_geometry, local_proj_stack,
+                                     local_volume);
     });
 }
 

examples/dose_deposition.cpp

@@ -6,7 +6,7 @@ int main() {
     int k = 256;
     auto v = tomo::volume<2_D, T>(k);
     auto f = tomo::modified_shepp_logan_phantom<T>(v);
-    auto g = tomo::geometry::parallel<2_D, T>(v, k, k);
+    auto g = tomo::geometry::parallel<2_D, T>(v, k);
 
     auto proj = tomo::dim::linear<2_D, T>(v);
     auto dose = tomo::image<2_D, T>(v);

examples/geometry_benchmark.cpp

@@ -106,7 +106,7 @@ int main(int argc, char* argv[]) {
     std::vector<std::thread> threads;
     // 'gather results'
     threads.emplace_back(partition_test<D, T>, "parallel",
-                         tomo::geometry::parallel<3_D, T>(v, k, k),
+                         tomo::geometry::parallel<3_D, T>(v, k),
                          std::ref(table), v, p, e);
 
     // threads.emplace_back(

examples/minimal_example.cpp

@@ -1,12 +1,13 @@
 #include "tomos/tomos.hpp"
+using namespace tomo;
 
 int main() {
     using T = float;
-    constexpr tomo::dimension D = 2_D;
+    constexpr dimension D = 2_D;
 
     int size = 128;
     auto v = tomo::volume<D, T>(size);
-    auto g = tomo::geometry::parallel<D, T>(v, size, size);
+    auto g = tomo::geometry::parallel<D, T>(v, size);
     auto f = tomo::modified_shepp_logan_phantom<T>(v);
     auto k = tomo::dim::joseph<D, T>(v);
     auto p = tomo::forward_projection<D, T>(f, g, k);

examples/read_geometry_test.cpp

@@ -1,6 +1,5 @@
 #include "tomos/tomos.hpp"
 #include "tomos/util/plotter.hpp"
-#include "tomos/util/rescale.hpp"
 
 #include "fmt/format.h"
 
@@ -34,18 +33,16 @@ int main(int argc, char* argv[]) {
 
         // std::cout << "done loading, now should downscale\n";
         proj_plotter.plot(problem.projection_stack.get_projection(0));
-        auto scaled_problem =
-            tomo::rescale<D, T>(problem, {64, 64}, {64, 64, 64}, 64);
 
-        auto proj = tomo::dim::joseph<D, T>(scaled_problem.object_volume);
+        auto proj = tomo::dim::joseph<D, T>(problem.object_volume);
 
         fmt::print("STARTING SIRT\n");
         auto z = tomo::reconstruction::sirt(
-            scaled_problem.object_volume, *scaled_problem.acquisition_geometry,
-            proj, scaled_problem.projection_stack, opt.beta, opt.iterations, {[&](tomo::image<D, T>& image) { plotter.plot(image); }});
+            problem.object_volume, *problem.acquisition_geometry,
+            proj, problem.projection_stack, opt.beta, opt.iterations, {[&](tomo::image<D, T>& image) { plotter.plot(image); }});
         fmt::print("SIRT\n");
         tomo::ascii_plot(z);
-        proj_plotter.plot(scaled_problem.projection_stack.get_projection(0));
+        proj_plotter.plot(problem.projection_stack.get_projection(0));
         plotter.plot(z);
     } catch (const std::exception& e) {
         std::cout << "Reading configuration failed: " << e.what() << "\n";

examples/reconstruction.cpp

@@ -10,7 +10,7 @@ using T = float;
 template <tomo::dimension D>
 void run(tomo::util::args opt) {
     auto v = tomo::volume<D, T>(opt.k);
-    auto g = tomo::geometry::parallel<D, T>(v, opt.k, opt.k);
+    auto g = tomo::geometry::parallel<D, T>(v, opt.k);
     auto f = tomo::modified_shepp_logan_phantom<T>(v);
 
     tomo::ascii_plot(f);

include/tomos/algorithms/sart.hpp

@@ -33,7 +33,7 @@ image<D, T> sart(const volume<D, T>& v, const tomo::geometry::base<D, T>& g,
     image<D, T> f(v);
 
     // the size of a single block
-    int k = g.detector_pixel_count();
+    int k = math::reduce<D -1>(g.projection_shape(0));
 
     // compute $w_i \cdot w_i$
     std::vector<T> w_norms(g.lines());

include/tomos/geometries/cone.hpp

@@ -68,17 +68,6 @@ class cone_beam : public trajectory<3_D, T> {
                 apply_rotation_(detector_tilt_[1], projection)};
     }
 
-    void set_projections(int projection_count) override {
-        angles_.clear();
-
-        T delta = ((T)2.0 * math::pi<T>) / projection_count;
-        for (int i = 0; i < projection_count; ++i) {
-            angles_.push_back(delta * i);
-        }
-
-        trajectory<3_D, T>::set_projections(projection_count);
-    }
-
   private:
     math::vec<3_D, T> source_position_;
     math::vec<3_D, T> detector_position_;
@@ -94,7 +83,6 @@ class cone_beam : public trajectory<3_D, T> {
 
     inline math::vec<3_D, T> apply_rotation_(math::vec<3_D, T> location,
                                              int projection) const {
-        // here we rotate around z
         return math::rotate(location, math::standard_basis<3_D, T>(2),
                             angles_[projection]);
     }

include/tomos/geometries/helical_cone_beam.hpp

@@ -73,17 +73,6 @@ class helical_cone_beam : public trajectory<3_D, T> {
                 apply_rotation_(detector_tilt_[1], projection)};
     }
 
-    void set_projections(int projection_count) override {
-        angles_.clear();
-
-        T delta = ((T)2.0 * math::pi<T>) / projection_count;
-        for (int i = 0; i < projection_count; ++i) {
-            angles_.push_back(delta * i);
-        }
-
-        trajectory<3_D, T>::set_projections(projection_count);
-    }
-
   private:
     math::vec<3_D, T> source_position_;
     math::vec<3_D, T> detector_position_;