diff --git a/src/aux.hpp b/src/aux.hpp index a815419..f720dc1 100644 --- a/src/aux.hpp +++ b/src/aux.hpp @@ -190,14 +190,13 @@ std::tuple cross_section_slow( * connected component, so pre-filtering must be performed to * ensure a match. */ -float cross_sectional_area_slow( +std::tuple cross_sectional_area_slow( const uint8_t* binimg, const uint64_t sx, const uint64_t sy, const uint64_t sz, const float px, const float py, const float pz, const float nx, const float ny, const float nz, - const float wx, const float wy, const float wz, - uint8_t &contact = _dummy_contact + const float wx, const float wy, const float wz ) { const Vec3 pos(px, py, pz); @@ -208,7 +207,7 @@ float cross_sectional_area_slow( || rpos.y < 0 || rpos.y >= sy || rpos.z < 0 || rpos.z >= sz ) { - return 0.0; + return std::make_tuple(0.0, 0); } const Vec3 anisotropy(wx, wy, wz); @@ -233,7 +232,7 @@ float cross_sectional_area_slow( : 1.0 / projections[i]; } - contact = 0; + uint8_t contact = 0; for (uint64_t z = 0; z < sz; z++) { for (uint64_t y = 0; y < sy; y++) { @@ -263,7 +262,7 @@ float cross_sectional_area_slow( } } - return area; + return std::make_tuple(area, contact); } }; diff --git a/src/fastxs3d.cpp b/src/fastxs3d.cpp index bc0fe13..7a9713d 100644 --- a/src/fastxs3d.cpp +++ b/src/fastxs3d.cpp @@ -84,7 +84,8 @@ auto calculate_area( const py::array_t &point, const py::array_t &normal, const py::array_t &anisotropy, - const bool slow_method + const bool slow_method, + const bool use_persistent_data ) { const uint64_t sx = binimg.shape()[0]; const uint64_t sy = binimg.ndim() < 2 @@ -94,31 +95,25 @@ auto calculate_area( ? 1 : binimg.shape()[2]; - uint8_t contact = false; - float area = 0; - if (slow_method) { - area = xs3d::cross_sectional_area_slow( + return xs3d::cross_sectional_area_slow( binimg.data(), sx, sy, sz, point.at(0), point.at(1), point.at(2), normal.at(0), normal.at(1), normal.at(2), - anisotropy.at(0), anisotropy.at(1), anisotropy.at(2), - contact + anisotropy.at(0), anisotropy.at(1), anisotropy.at(2) ); } else { - area = xs3d::cross_sectional_area( + return xs3d::cross_sectional_area( binimg.data(), sx, sy, sz, point.at(0), point.at(1), point.at(2), normal.at(0), normal.at(1), normal.at(2), anisotropy.at(0), anisotropy.at(1), anisotropy.at(2), - contact + use_persistent_data ); } - - return std::tuple(area, contact); } auto projection( @@ -225,4 +220,7 @@ PYBIND11_MODULE(fastxs3d, m) { m.def("projection", &projection, "Project a cross section of a 3D image onto a 2D plane"); m.def("section", §ion, "Return a floating point image that shows the voxels contributing area to a cross section."); m.def("area", &calculate_area, "Find the cross sectional area for a given binary image, point, and normal vector."); + + m.def("set_shape", &xs3d::set_shape, "Accelerate the area function across many evaluation points by saving some attributes of the input shape upfront. Call clear_shape when you are done."); + m.def("clear_shape", &xs3d::clear_shape, "Delete the data that was persisted by set_shape."); } diff --git a/src/xs3d.hpp b/src/xs3d.hpp index 807eab4..b52401d 100644 --- a/src/xs3d.hpp +++ b/src/xs3d.hpp @@ -20,6 +20,83 @@ namespace { static uint8_t _dummy_contact = false; +// half rounded up +uint64_t _h(uint64_t a) { + return ((a+1) >> 1); +}; + +uint8_t compute_cube( + const uint8_t* binimg, + const uint64_t sx, const uint64_t sy, const uint64_t sz, + const uint64_t x, const uint64_t y, const uint64_t z +) { + const uint64_t sxy = sx * sy; + const uint64_t loc = x + sx * (y + sy * z); + + const uint64_t x_valid = (x < sx - 1); + const uint64_t y_valid = (y < sy - 1); + const uint64_t z_valid = (z < sz - 1); + + return static_cast( + (binimg[loc] > 0) + | ((x_valid && (binimg[loc+1] > 0)) << 1) + | ((y_valid && (binimg[loc+sx] > 0)) << 2) + | (((x_valid && y_valid) && (binimg[loc+sx+1] > 0)) << 3) + | ((z_valid && (binimg[loc+sxy] > 0)) << 4) + | (((x_valid && z_valid) && (binimg[loc+sxy+1] > 0)) << 5) + | (((y_valid && z_valid) && (binimg[loc+sxy+sx] > 0)) << 6) + | (((x_valid && y_valid && z_valid) && (binimg[loc+sxy+sx+1] > 0)) << 7) + ); +} + +struct PersistentShapeManager { + std::vector visited; + uint64_t sx, sy, sz; + uint8_t color; + + PersistentShapeManager() { + sx = 0; sy = 0; sz = 0; + color = 0; + } + + void init(const uint64_t _sx, const uint64_t _sy, const uint64_t _sz) { + sx = _sx; + sy = _sy; + sz = _sz; + visited.resize(this->eighth_voxels()); + } + + void maybe_grow(const uint64_t _sx, const uint64_t _sy, const uint64_t _sz) { + if (sx * sy * sz < _sx * _sy * _sz) { + init(_sx, _sy, _sz); + } + } + + void clear() { + sx = 0; + sy = 0; + sz = 0; + color = 0; + visited = std::vector(); + } + + uint64_t eighth_voxels() { + return _h(sx) * _h(sy) * _h(sz); + } + + uint8_t next_color() { + color++; + if (color == 0) { + std::fill(visited.begin(), visited.end(), 0); + color = 1; + } + + return color; + } +}; + +PersistentShapeManager persistent_data; + const Vec3 ihat = Vec3(1,0,0); const Vec3 jhat = Vec3(0,1,0); const Vec3 khat = Vec3(0,0,1); @@ -459,26 +536,6 @@ float calc_area_at_point( return subtotal; } -uint8_t compute_cube( - const uint8_t* binimg, - const uint64_t sx, const uint64_t sy, const uint64_t sz, - const uint64_t x, const uint64_t y, const uint64_t z -) { - const uint64_t sxy = sx * sy; - const uint64_t loc = x + sx * (y + sy * z); - - return static_cast( - (binimg[loc] > 0) - | (((x < sx - 1) && (binimg[loc+1] > 0)) << 1) - | (((y < sy - 1) && (binimg[loc+sx] > 0)) << 2) - | (((x < sx - 1 && y < sy - 1) && (binimg[loc+sx+1] > 0)) << 3) - | (((z < sz - 1) && (binimg[loc+sxy] > 0)) << 4) - | (((x < sx - 1 && z < sz - 1) && (binimg[loc+sxy+1] > 0)) << 5) - | (((y < sy - 1 && z < sz - 1) && (binimg[loc+sxy+sx] > 0)) << 6) - | (((x < sx - 1 && y < sy - 1 && z < sz - 1) && (binimg[loc+sxy+sx+1] > 0)) << 7) - ); -} - bool is_26_connected( const uint8_t center, const uint8_t candidate, const int x, const int y, const int z @@ -642,8 +699,8 @@ float robust_calc_area_at_point_2x2x2( static_cast(delta.y) + sy * static_cast(delta.z) ); - const uint64_t ccl_loc = (static_cast(delta.x) >> 1) + ((sx+1) >> 1) * ( - (static_cast(delta.y) >> 1) + ((sy+1) >> 1) * (static_cast(delta.z) >> 1) + const uint64_t ccl_loc = (static_cast(delta.x) >> 1) + _h(sx) * ( + (static_cast(delta.y) >> 1) + _h(sy) * (static_cast(delta.z) >> 1) ); if (!binimg[loc] || ccl[ccl_loc]) { @@ -670,17 +727,99 @@ float robust_calc_area_at_point_2x2x2( return subtotal; } -float cross_sectional_area_helper_2x2x2( +float robust_calc_area_at_point_2x2x2_persistent_data( + const uint8_t* binimg, + const uint64_t sx, const uint64_t sy, const uint64_t sz, + const Vec3& cur, const Vec3& pos, + const Vec3& normal, const Vec3& anisotropy, + std::vector& pts, + const std::vector& projections, + const std::vector& inv_projections +) { + + uint64_t x = static_cast(cur.x) & ~static_cast(1); + uint64_t y = static_cast(cur.y) & ~static_cast(1); + uint64_t z = static_cast(cur.z) & ~static_cast(1); + + float subtotal = 0.0; + + float xs = (cur.x - 2) >= 0 ? -2 : 0; + float ys = (cur.y - 2) >= 0 ? -2 : 0; + float zs = (cur.z - 2) >= 0 ? -2 : 0; + + float xe = (cur.x + 2) < sx ? 2 : 0; + float ye = (cur.y + 2) < sy ? 2 : 0; + float ze = (cur.z + 2) < sz ? 2 : 0; + + // only need to check around the current voxel if + // there's a possibility that there is a gap due + // to basis vector motion. If the normal is axis + // aligned to x, y, or z, there will be no gap. + if (normal.is_axis_aligned()) { + xs = 0; + ys = 0; + zs = 0; + + xe = 0; + ye = 0; + ze = 0; + } + + const uint8_t center = compute_cube(binimg, sx, sy, sz, x, y, z); + std::vector& visited = persistent_data.visited; + const uint8_t color = persistent_data.color; + + for (int64_t zi = zs; zi <= ze; zi += 2) { + for (int64_t yi = ys; yi <= ye; yi += 2) { + for (int64_t xi = xs; xi <= xe; xi += 2) { + + Vec3 delta(xi,yi,zi); + delta += cur; + + const uint64_t loc = static_cast(delta.x) + sx * ( + static_cast(delta.y) + sy * static_cast(delta.z) + ); + + const uint64_t visited_loc = (static_cast(delta.x) >> 1) + _h(sx) * ( + (static_cast(delta.y) >> 1) + _h(sy) * (static_cast(delta.z) >> 1) + ); + + if (!binimg[loc] || visited[visited_loc] == color) { + continue; + } + + visited[visited_loc] = color; + + const uint8_t candidate = compute_cube(binimg, sx, sy, sz, x + xi, y + yi, z + zi); + + if (is_26_connected(center, candidate, xi, yi, zi)) { + subtotal += calc_area_at_point_2x2x2( + candidate, + sx, sy, sz, + delta, pos, normal, anisotropy, + pts, + projections, inv_projections + ); + } + } + } + } + + return subtotal; +} + +std::tuple cross_sectional_area_helper_2x2x2( const uint8_t* binimg, const uint64_t sx, const uint64_t sy, const uint64_t sz, const Vec3& pos, // plane position const Vec3& normal, // plane normal vector - const Vec3& anisotropy, // anisotropy - uint8_t& contact + const Vec3& anisotropy ) { const uint64_t grid_size = std::max(((sx+1)>>1) * ((sy+1)>>1) * ((sz+1)>>1), static_cast(1)); std::vector ccl(grid_size); + uint8_t contact = 0; + // rational approximation of sqrt(3) is 97/56 // more reliable behavior across compilers/architectures uint64_t plane_size = 2 * 97 * std::max(std::max(sx,sy), sz) / 56 + 1; @@ -820,9 +959,164 @@ float cross_sectional_area_helper_2x2x2( ); } - return total; + return std::make_tuple(total, contact); +} + +std::tuple cross_sectional_area_helper_2x2x2_persistent_data( + const uint8_t* binimg, + const uint64_t sx, const uint64_t sy, const uint64_t sz, + const Vec3& pos, // plane position + const Vec3& normal, // plane normal vector + const Vec3& anisotropy +) { + persistent_data.maybe_grow(sx, sy, sz); + persistent_data.next_color(); + + uint8_t contact = 0; + + // rational approximation of sqrt(3) is 97/56 + // more reliable behavior across compilers/architectures + uint64_t plane_size = 2 * 97 * std::max(std::max(sx,sy), sz) / 56 + 1; + + // maximum possible size of plane + uint64_t psx = plane_size; + uint64_t psy = psx; + + std::vector visited(psx * psy); + + Vec3 basis1 = normal.cross(ihat); + if (basis1.is_null()) { + basis1 = normal.cross(jhat); + } + basis1 /= basis1.norm(); + + Vec3 basis2 = normal.cross(basis1); + basis2 /= basis2.norm(); + + uint64_t plane_pos_x = plane_size / 2; + uint64_t plane_pos_y = plane_size / 2; + + uint64_t ploc = plane_pos_x + psx * plane_pos_y; + + std::stack stack; + stack.push(ploc); + + float total = 0.0; + + std::vector pts; + pts.reserve(6); + + const std::vector projections = { + ihat.dot(normal), + jhat.dot(normal), + khat.dot(normal) + }; + + std::vector inv_projections(3); + for (int i = 0; i < 3; i++) { + inv_projections[i] = (projections[i] == 0) + ? 0 + : 1.0 / projections[i]; + } + + const float sxf = static_cast(sx) - 0.5; + const float syf = static_cast(sy) - 0.5; + const float szf = static_cast(sz) - 0.5; + + while (!stack.empty()) { + ploc = stack.top(); + stack.pop(); + + if (visited[ploc]) { + continue; + } + + visited[ploc] = true; + + uint64_t y = ploc / psx; + uint64_t x = ploc - y * psx; + + float dx = static_cast(x) - static_cast(plane_pos_x); + float dy = static_cast(y) - static_cast(plane_pos_y); + + Vec3 cur = pos + basis1 * dx + basis2 * dy; + + if (cur.x < -0.5 || cur.y < -0.5 || cur.z < -0.5) { + continue; + } + else if (cur.x >= sxf || cur.y >= syf || cur.z >= szf) { + continue; + } + + cur = cur.round(); + + uint64_t loc = ( + static_cast(cur.x) + + sx * ( + static_cast(cur.y) + + sy * static_cast(cur.z) + ) + ); + + if (!binimg[loc]) { + continue; + } + + contact |= (cur.x < 1); // -x + contact |= (cur.x >= sx - 1.5) << 1; // +x + contact |= (cur.y < 1) << 2; // -y + contact |= (cur.y >= sy - 1.5) << 3; // +y + contact |= (cur.z < 1) << 4; // -z + contact |= (cur.z >= sz - 1.5) << 5; // +z + + uint64_t up = ploc - psx; + uint64_t down = ploc + psx; + uint64_t left = ploc - 1; + uint64_t right = ploc + 1; + + uint64_t upleft = ploc - psx - 1; + uint64_t downleft = ploc + psx - 1; + uint64_t upright = ploc - psx + 1; + uint64_t downright = ploc + psx + 1; + + if (x > 0 && !visited[left]) { + stack.push(left); + } + if (x < psx - 1 && !visited[right]) { + stack.push(right); + } + if (y > 0 && !visited[up]) { + stack.push(up); + } + if (y < psy - 1 && !visited[down]) { + stack.push(down); + } + + if (x > 0 && y > 0 && !visited[upleft]) { + stack.push(upleft); + } + if (x < psx - 1 && y > 0 && !visited[upright]) { + stack.push(upright); + } + if (x > 0 && y < psy - 1 && !visited[downleft]) { + stack.push(downleft); + } + if (x < psx - 1 && y < psy - 1 && !visited[downright]) { + stack.push(downright); + } + + total += robust_calc_area_at_point_2x2x2_persistent_data( + binimg, + sx, sy, sz, + cur, pos, normal, anisotropy, + pts, projections, inv_projections + ); + } + + return std::make_tuple(total, contact); } + float cross_sectional_area_helper( const uint8_t* binimg, const uint64_t sx, const uint64_t sy, const uint64_t sz, @@ -1000,24 +1294,34 @@ struct Bbox2d { } }; -float cross_sectional_area( +void set_shape( + const uint64_t sx, const uint64_t sy, const uint64_t sz +) { + persistent_data.init(sx, sy, sz); +} + +void clear_shape() { + persistent_data.clear(); +} + +std::tuple cross_sectional_area( const uint8_t* binimg, const uint64_t sx, const uint64_t sy, const uint64_t sz, const float px, const float py, const float pz, const float nx, const float ny, const float nz, const float wx, const float wy, const float wz, - uint8_t &contact = _dummy_contact + const bool use_persistent_data = false ) { if (px < 0 || px >= sx) { - return 0.0; + return std::make_tuple(0.0, 0); } else if (py < 0 || py >= sy) { - return 0.0; + return std::make_tuple(0.0, 0); } else if (pz < 0 || pz >= sz) { - return 0.0; + return std::make_tuple(0.0, 0); } const Vec3 pos(px, py, pz); @@ -1028,7 +1332,7 @@ float cross_sectional_area( || rpos.y < 0 || rpos.y >= sy || rpos.z < 0 || rpos.z >= sz ) { - return 0.0; + return std::make_tuple(0.0, 0); } uint64_t loc = static_cast(rpos.x) + sx * ( @@ -1036,22 +1340,30 @@ float cross_sectional_area( ); if (loc < 0 || loc >= sx * sy * sz) { - return 0.0; + return std::make_tuple(0.0, 0); } else if (!binimg[loc]) { - return 0.0; + return std::make_tuple(0.0, 0); } const Vec3 anisotropy(wx, wy, wz); Vec3 normal(nx, ny, nz); normal /= normal.norm(); - return cross_sectional_area_helper_2x2x2( - binimg, - sx, sy, sz, - pos, normal, anisotropy, - contact - ); + if (use_persistent_data) { + return cross_sectional_area_helper_2x2x2_persistent_data( + binimg, + sx, sy, sz, + pos, normal, anisotropy + ); + } + else { + return cross_sectional_area_helper_2x2x2( + binimg, + sx, sy, sz, + pos, normal, anisotropy + ); + } } std::tuple cross_section( diff --git a/xs3d/__init__.py b/xs3d/__init__.py index 48ed874..e9775ad 100644 --- a/xs3d/__init__.py +++ b/xs3d/__init__.py @@ -15,6 +15,7 @@ def cross_sectional_area( anisotropy:Optional[VECTOR_T] = None, return_contact:bool = False, slow_method:bool = False, + use_persistent_data:bool = False, ) -> Union[float, tuple[float, int]]: """ Find the cross sectional area for a given binary image, @@ -47,6 +48,12 @@ def cross_sectional_area( all locations are visited. Does not restrict analysis to a single connected component. + use_persistent_data: Use a pre-allocated buffer for + internally tacking visited positions. You allocate the + buffer with xs3d.set_shape and clear it with xs3d.clear_shape. + This can save about 20% of the time when repeatedly + analyzing a shape. + Returns: physical area covered by the section plane """ if anisotropy is None: @@ -75,7 +82,7 @@ def cross_sectional_area( area, contact = fastxs3d.area( binimg.view(np.uint8), pos, normal, anisotropy, - slow_method + slow_method, use_persistent_data, ) else: raise ValueError("dimensions not supported") @@ -224,7 +231,20 @@ def slice( crop ) +def set_shape(image:npt.NDArray[np.integer]): + """ + Allocate a buffer appropriately sized to this image for internal use. + This can accelerate the area calculation if there are repeated queries + against the same shape. + + It is very important that the persisted shape match the current input + or memory issues can result. + """ + fastxs3d.set_shape(image.shape[0], image.shape[1], image.shape[2]) +def clear_shape(): + """Free the persisted memory from set_shape.""" + fastxs3d.clear_shape()