Skip to content

feat: support for multilabel images #18

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
william-silversmith merged 4 commits into from
Oct 17, 2025
Merged

feat: support for multilabel images #18

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
0114518
feat: multilabel support
william-silversmith Oct 17, 2025
b394ae5
fix: boolean images are no longer required
william-silversmith Oct 17, 2025
08badff
fix: set setid to 1 only for binary images
william-silversmith Oct 17, 2025
3bc493b
docs: update to note that 1 for binary images
william-silversmith Oct 17, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions src/auxiliary.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -190,8 +190,9 @@ std::tuple<float*, uint8_t> cross_section_slow(
* connected component, so pre-filtering must be performed to
* ensure a match.
*/
template <typename LABEL>
std::tuple<float, uint8_t> cross_sectional_area_slow(
const uint8_t* binimg,
const LABEL* labels, const LABEL segid,
const uint64_t sx, const uint64_t sy, const uint64_t sz,

const float px, const float py, const float pz,
Expand Down Expand Up @@ -239,7 +240,7 @@ std::tuple<float, uint8_t> cross_sectional_area_slow(
for (uint64_t x = 0; x < sx; x++) {
uint64_t loc = x + sx * (y + sy * z);

if (!binimg[loc]) {
if (labels[loc] != segid) {
continue;
}

Expand Down
41 changes: 30 additions & 11 deletions src/fastxs3d.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -79,21 +79,23 @@ auto section(
return std::make_tuple(arr, std::get<1>(tup));
}

template <typename LABEL>
auto calculate_area(
const py::array_t<uint8_t> &binimg,
const py::array_t<LABEL> &image,
const LABEL segid,
const py::array_t<float> &point,
const py::array_t<float> &normal,
const py::array_t<float> &anisotropy,
const bool slow_method,
const bool use_persistent_data
) {
const uint64_t sx = binimg.shape()[0];
const uint64_t sy = binimg.ndim() < 2
const uint64_t sx = image.shape()[0];
const uint64_t sy = image.ndim() < 2
? 1
: binimg.shape()[1];
const uint64_t sz = binimg.ndim() < 3
: image.shape()[1];
const uint64_t sz = image.ndim() < 3
? 1
: binimg.shape()[2];
: image.shape()[2];

if (point.size() < 3) {
throw py::value_error("point array must have at least 3 elements");
Expand All @@ -112,17 +114,17 @@ auto calculate_area(
}

if (slow_method) {
return xs3d::cross_sectional_area_slow(
binimg.data(),
return xs3d::cross_sectional_area_slow<LABEL>(
image.data(), segid,
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)
);
}
else {
return xs3d::cross_sectional_area(
binimg.data(),
return xs3d::cross_sectional_area<LABEL>(
image.data(), segid,
sx, sy, sz,
point.at(0), point.at(1), point.at(2),
normal.at(0), normal.at(1), normal.at(2),
Expand Down Expand Up @@ -231,12 +233,29 @@ auto projection(
}
}


#define REGISTER_CALCULATE_AREA_FOR_TYPE(T) \
m.def("area", &calculate_area<T>, "Find the cross sectional area for a given image, point, and normal vector.");

PYBIND11_MODULE(fastxs3d, m) {
m.doc() = "Finding cross sectional area in 3D voxelized images.";
m.def("projection", &projection, "Project a cross section of a 3D image onto a 2D plane");
m.def("section", &section, "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.");

REGISTER_CALCULATE_AREA_FOR_TYPE(bool)
REGISTER_CALCULATE_AREA_FOR_TYPE(uint8_t)
REGISTER_CALCULATE_AREA_FOR_TYPE(uint16_t)
REGISTER_CALCULATE_AREA_FOR_TYPE(uint32_t)
REGISTER_CALCULATE_AREA_FOR_TYPE(uint64_t)
REGISTER_CALCULATE_AREA_FOR_TYPE(int8_t)
REGISTER_CALCULATE_AREA_FOR_TYPE(int16_t)
REGISTER_CALCULATE_AREA_FOR_TYPE(int32_t)
REGISTER_CALCULATE_AREA_FOR_TYPE(int64_t)
REGISTER_CALCULATE_AREA_FOR_TYPE(float)
REGISTER_CALCULATE_AREA_FOR_TYPE(double)

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.");
}

#undef REGISTER_CALCULATE_AREA_FOR_TYPE
81 changes: 44 additions & 37 deletions src/xs3d.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,8 +23,9 @@ uint64_t _h(uint64_t a) {
return ((a+1) >> 1);
};

template <typename LABEL>
uint8_t compute_cube(
const uint8_t* binimg,
const LABEL* labels, const LABEL segid,
const uint64_t sx, const uint64_t sy, const uint64_t sz,
const uint64_t x, const uint64_t y, const uint64_t z
) {
Expand All @@ -36,14 +37,14 @@ uint8_t compute_cube(
const uint64_t z_valid = (z < sz - 1);

return static_cast<uint8_t>(
(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)
(labels[loc] == segid)
| ((x_valid && (labels[loc+1] == segid)) << 1)
| ((y_valid && (labels[loc+sx] == segid)) << 2)
| (((x_valid && y_valid) && (labels[loc+sx+1] == segid)) << 3)
| ((z_valid && (labels[loc+sxy] == segid)) << 4)
| (((x_valid && z_valid) && (labels[loc+sxy+1] == segid)) << 5)
| (((y_valid && z_valid) && (labels[loc+sxy+sx] == segid)) << 6)
| (((x_valid && y_valid && z_valid) && (labels[loc+sxy+sx+1] == segid)) << 7)
);
}

Expand Down Expand Up @@ -439,8 +440,9 @@ float calc_area_at_point_2x2x2(
return area;
}

template <typename LABEL>
float calc_area_at_point(
const uint8_t* binimg,
const LABEL* labels, const LABEL segid,
std::vector<bool>& ccl,
const uint64_t sx, const uint64_t sy, const uint64_t sz,
const Vec3& cur, const Vec3& pos,
Expand Down Expand Up @@ -504,7 +506,7 @@ float calc_area_at_point(
if (loc < 0 || loc >= voxels) {
continue;
}
else if (!binimg[loc]) {
else if (labels[loc] != segid) {
continue;
}

Expand Down Expand Up @@ -645,8 +647,9 @@ bool is_26_connected(
}
}

template <typename LABEL>
float robust_calc_area_at_point_2x2x2(
const uint8_t* binimg,
const LABEL* labels, const LABEL segid,
std::vector<bool>& ccl,
const uint64_t sx, const uint64_t sy, const uint64_t sz,
const Vec3& cur, const Vec3& pos,
Expand Down Expand Up @@ -684,7 +687,7 @@ float robust_calc_area_at_point_2x2x2(
ze = 0;
}

const uint8_t center = compute_cube(binimg, sx, sy, sz, x, y, z);
const uint8_t center = compute_cube(labels, segid, sx, sy, sz, x, y, z);

for (int64_t zi = zs; zi <= ze; zi += 2) {
for (int64_t yi = ys; yi <= ye; yi += 2) {
Expand All @@ -701,13 +704,13 @@ float robust_calc_area_at_point_2x2x2(
(static_cast<uint64_t>(delta.y) >> 1) + _h(sy) * (static_cast<uint64_t>(delta.z) >> 1)
);

if (!binimg[loc] || ccl[ccl_loc]) {
if ((labels[loc] != segid) || ccl[ccl_loc]) {
continue;
}

ccl[ccl_loc] = true;

uint8_t candidate = compute_cube(binimg, sx, sy, sz, x + xi, y + yi, z + zi);
uint8_t candidate = compute_cube(labels, segid, sx, sy, sz, x + xi, y + yi, z + zi);

if (is_26_connected(center, candidate, xi, yi, zi)) {
subtotal += calc_area_at_point_2x2x2(
Expand All @@ -725,8 +728,9 @@ float robust_calc_area_at_point_2x2x2(
return subtotal;
}

template <typename LABEL>
float robust_calc_area_at_point_2x2x2_persistent_data(
const uint8_t* binimg,
const LABEL* labels, const LABEL segid,
const uint64_t sx, const uint64_t sy, const uint64_t sz,
const Vec3& cur, const Vec3& pos,
const Vec3& normal, const Vec3& anisotropy,
Expand Down Expand Up @@ -763,7 +767,7 @@ float robust_calc_area_at_point_2x2x2_persistent_data(
ze = 0;
}

const uint8_t center = compute_cube(binimg, sx, sy, sz, x, y, z);
const uint8_t center = compute_cube(labels, segid, sx, sy, sz, x, y, z);
std::vector<uint8_t>& visited = persistent_data.visited;
const uint8_t color = persistent_data.color;

Expand All @@ -782,13 +786,13 @@ float robust_calc_area_at_point_2x2x2_persistent_data(
(static_cast<uint64_t>(delta.y) >> 1) + _h(sy) * (static_cast<uint64_t>(delta.z) >> 1)
);

if (!binimg[loc] || visited[visited_loc] == color) {
if ((labels[loc] != segid) || 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);
const uint8_t candidate = compute_cube(labels, segid, sx, sy, sz, x + xi, y + yi, z + zi);

if (is_26_connected(center, candidate, xi, yi, zi)) {
subtotal += calc_area_at_point_2x2x2(
Expand All @@ -806,8 +810,9 @@ float robust_calc_area_at_point_2x2x2_persistent_data(
return subtotal;
}

template <typename LABEL>
std::tuple<float, uint8_t> cross_sectional_area_helper_2x2x2(
const uint8_t* binimg,
const LABEL* labels, const LABEL segid,
const uint64_t sx, const uint64_t sy, const uint64_t sz,
const Vec3& pos, // plane position
const Vec3& normal, // plane normal vector
Expand Down Expand Up @@ -902,7 +907,7 @@ std::tuple<float, uint8_t> cross_sectional_area_helper_2x2x2(
)
);

if (!binimg[loc]) {
if (labels[loc] != segid) {
continue;
}

Expand Down Expand Up @@ -949,8 +954,8 @@ std::tuple<float, uint8_t> cross_sectional_area_helper_2x2x2(
stack.push(downright);
}

total += robust_calc_area_at_point_2x2x2(
binimg, ccl,
total += robust_calc_area_at_point_2x2x2<LABEL>(
labels, segid, ccl,
sx, sy, sz,
cur, pos, normal, anisotropy,
pts, projections, inv_projections
Expand All @@ -960,8 +965,9 @@ std::tuple<float, uint8_t> cross_sectional_area_helper_2x2x2(
return std::make_tuple(total, contact);
}

template <typename LABEL>
std::tuple<float, uint8_t> cross_sectional_area_helper_2x2x2_persistent_data(
const uint8_t* binimg,
const LABEL* labels, const LABEL segid,
const uint64_t sx, const uint64_t sy, const uint64_t sz,
const Vec3& pos, // plane position
const Vec3& normal, // plane normal vector
Expand Down Expand Up @@ -1056,7 +1062,7 @@ std::tuple<float, uint8_t> cross_sectional_area_helper_2x2x2_persistent_data(
)
);

if (!binimg[loc]) {
if (labels[loc] != segid) {
continue;
}

Expand Down Expand Up @@ -1104,7 +1110,7 @@ std::tuple<float, uint8_t> cross_sectional_area_helper_2x2x2_persistent_data(
}

total += robust_calc_area_at_point_2x2x2_persistent_data(
binimg,
labels, segid,
sx, sy, sz,
cur, pos, normal, anisotropy,
pts, projections, inv_projections
Expand All @@ -1114,9 +1120,9 @@ std::tuple<float, uint8_t> cross_sectional_area_helper_2x2x2_persistent_data(
return std::make_tuple(total, contact);
}


template <typename LABEL>
float cross_sectional_area_helper(
const uint8_t* binimg,
const LABEL* labels, const LABEL segid,
const uint64_t sx, const uint64_t sy, const uint64_t sz,
const Vec3& pos, // plane position
const Vec3& normal, // plane normal vector
Expand Down Expand Up @@ -1208,7 +1214,7 @@ float cross_sectional_area_helper(
)
);

if (!binimg[loc]) {
if (labels[loc] != segid) {
continue;
}

Expand Down Expand Up @@ -1256,7 +1262,7 @@ float cross_sectional_area_helper(
}

total += calc_area_at_point(
binimg, ccl,
labels, segid, ccl,
sx, sy, sz,
cur, pos, normal, anisotropy,
pts, projections, inv_projections,
Expand Down Expand Up @@ -1302,8 +1308,9 @@ void clear_shape() {
persistent_data.clear();
}

template <typename LABEL>
std::tuple<float, uint8_t> cross_sectional_area(
const uint8_t* binimg,
const LABEL* labels, const LABEL segid,
const uint64_t sx, const uint64_t sy, const uint64_t sz,

const float px, const float py, const float pz,
Expand Down Expand Up @@ -1340,7 +1347,7 @@ std::tuple<float, uint8_t> cross_sectional_area(
if (loc < 0 || loc >= sx * sy * sz) {
return std::make_tuple(0.0, 0);
}
else if (!binimg[loc]) {
else if (labels[loc] != segid) {
return std::make_tuple(0.0, 0);
}

Expand All @@ -1349,15 +1356,15 @@ std::tuple<float, uint8_t> cross_sectional_area(
normal /= normal.norm();

if (use_persistent_data) {
return cross_sectional_area_helper_2x2x2_persistent_data(
binimg,
return cross_sectional_area_helper_2x2x2_persistent_data<LABEL>(
labels, segid,
sx, sy, sz,
pos, normal, anisotropy
);
}
else {
return cross_sectional_area_helper_2x2x2(
binimg,
return cross_sectional_area_helper_2x2x2<LABEL>(
labels, segid,
sx, sy, sz,
pos, normal, anisotropy
);
Expand Down Expand Up @@ -1416,7 +1423,7 @@ std::tuple<float*, uint8_t> cross_section(
normal /= normal.norm();

cross_sectional_area_helper(
binimg,
binimg, /*segid=*/static_cast<uint8_t>(1),
sx, sy, sz,
pos, normal, anisotropy,
contact, plane_visualization
Expand Down
Loading