diff --git a/src/fastxs3d.cpp b/src/fastxs3d.cpp index 59165f5..7e8416d 100644 --- a/src/fastxs3d.cpp +++ b/src/fastxs3d.cpp @@ -83,40 +83,20 @@ auto projection( ? 1 : labels.shape()[2]; - float wx = anisotropy.at(0); - float wy = anisotropy.at(1); - float wz = anisotropy.at(2); - - float minval = std::min(std::min(wx,wy), wz); - wx /= minval; wy /= minval; wz /= minval; - float maxval = std::max(std::max(std::abs(wx), std::abs(wy)), std::abs(wz)); - - const uint64_t distortion = static_cast(ceil(maxval)); - - // rational approximation of sqrt(3) is 97/56 - // result is more likely to be same across compilers - uint64_t psx = distortion * 2 * 97 * std::max(std::max(sx,sy), sz) / 56 + 1; - uint64_t pvoxels = psx * psx; - - py::array arr; - auto projectionfn = [&](auto dtype) { - arr = py::array_t({ psx, psx }); - auto out = reinterpret_cast(arr.request().ptr); auto data = reinterpret_cast(labels.request().ptr); - std::fill(out, out + pvoxels, 0); - - std::tuple tup = xs3d::cross_section_projection( + auto tup = xs3d::cross_section_projection( 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), - standardize_basis, - out + standardize_basis ); - xs3d::Bbox2d bbox = std::get<1>(tup); + auto out = std::get<0>(tup); + uint64_t psx = std::get<1>(tup); + xs3d::Bbox2d bbox = std::get<2>(tup); bbox.x_max++; bbox.y_max++; @@ -133,6 +113,8 @@ auto projection( } } + delete[] out; + return cutout.view(py::str(labels.dtype())); }; diff --git a/src/xs3d.hpp b/src/xs3d.hpp index 70f9a05..d78fcdc 100644 --- a/src/xs3d.hpp +++ b/src/xs3d.hpp @@ -10,10 +10,8 @@ #include #include -namespace { - -static uint8_t _dummy_contact = false; +namespace xs3d { class Vec3 { public: float x, y, z; @@ -137,6 +135,35 @@ class Vec3 { } }; +struct Bbox2d { + int64_t x_min, x_max; + int64_t y_min, y_max; + Bbox2d() : x_min(0), x_max(0), y_min(0), y_max(0) {}; + Bbox2d(int64_t x_min, int64_t x_max, int64_t y_min, int64_t y_max) + : x_min(x_min), x_max(x_max), y_min(y_min), y_max(y_max) {}; + + int64_t sx() const { + return x_max - x_min; + } + int64_t sy() const { + return y_max - y_min; + } + int64_t pixels() const { + return sx() * sy(); + } + void print() const { + printf("Bbox2d(%lld, %lld, %lld, %lld)\n", x_min, x_max, y_min, y_max); + } +}; + +}; + +namespace { + +using namespace xs3d; + +static uint8_t _dummy_contact = false; + const Vec3 ihat = Vec3(1,0,0); const Vec3 jhat = Vec3(0,1,0); const Vec3 khat = Vec3(0,0,1); @@ -603,27 +630,6 @@ float cross_sectional_area_helper( namespace xs3d { -struct Bbox2d { - int64_t x_min, x_max; - int64_t y_min, y_max; - Bbox2d() : x_min(0), x_max(0), y_min(0), y_max(0) {}; - Bbox2d(int64_t x_min, int64_t x_max, int64_t y_min, int64_t y_max) - : x_min(x_min), x_max(x_max), y_min(y_min), y_max(y_max) {}; - - int64_t sx() const { - return x_max - x_min; - } - int64_t sy() const { - return y_max - y_min; - } - int64_t pixels() const { - return sx() * sy(); - } - void print() const { - printf("Bbox2d(%llu, %llu, %llu, %llu)\n", x_min, x_max, y_min, y_max); - } -}; - float cross_sectional_area( const uint8_t* binimg, const uint64_t sx, const uint64_t sy, const uint64_t sz, @@ -714,8 +720,13 @@ float* cross_section( } std::tuple create_orthonormal_basis( - const Vec3& normal, const bool positive_basis + const Vec3& normal, + const Vec3& anisotropy, + const bool positive_basis ) { + Vec3 distortion = anisotropy / anisotropy.min(); + distortion = Vec3(1,1,1) / distortion; + Vec3 basis1 = normal.cross(jhat); if (basis1.is_null()) { basis1 = normal.cross(ihat); @@ -755,128 +766,122 @@ std::tuple create_orthonormal_basis( } } - return std::tuple(basis1, basis2); + return std::tuple(basis1 * distortion, basis2 * distortion); +} + +Bbox2d compute_slice_plane( + const Vec3& pos, + const Vec3& basis1, const Vec3& basis2, + const uint64_t sx, const uint64_t sy, const uint64_t sz +) { + + auto minfn = [&](const Vec3& basis) { + float dxx = -pos.x * sqrt(3) / basis.x; + float dxy = -pos.y * sqrt(3) / basis.y; + float dxz = -pos.z * sqrt(3) / basis.z; + + dxx = std::isinf(dxx) ? INFINITY : dxx; + dxy = std::isinf(dxy) ? INFINITY : dxy; + dxz = std::isinf(dxz) ? INFINITY : dxz; + + return floor(std::min(std::min(dxx, dxy), dxz)); + }; + + auto maxfn = [&](const Vec3& basis) { + float dyx = (sqrt(3) * sx-pos.x) / basis.x; + float dyy = (sqrt(3) * sy-pos.y) / basis.y; + float dyz = (sqrt(3) * sz-pos.z) / basis.z; + + dyx = std::isinf(dyx)? -INFINITY : dyx; + dyy = std::isinf(dyy)? -INFINITY : dyy; + dyz = std::isinf(dyz)? -INFINITY : dyz; + + return ceil(std::max(std::max(dyx, dyy), dyz)); + }; + + float dx_min = minfn(basis1) - 1; + float dx_max = maxfn(basis1) + 1; + + float dy_min = minfn(basis2) - 1; + float dy_max = maxfn(basis2) + 1; + + return Bbox2d(dx_min, dx_max, dy_min, dy_max); } template -std::tuple cross_section_projection( +std::tuple cross_section_projection( const LABEL* labels, 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, - const bool positive_basis, - LABEL* out = NULL + const bool positive_basis ) { + const Vec3 anisotropy(wx, wy, wz); + const Vec3 pos(px, py, pz); - Vec3 anisotropy(wx, wy, wz); - anisotropy /= anisotropy.min(); - const uint64_t distortion = static_cast(ceil( - anisotropy.abs().max() - )); - anisotropy = Vec3(1,1,1) / anisotropy; - - // maximum possible size of plane - // rational approximation of sqrt(3) is 97/56 - const uint64_t psx = (distortion * 2 * 97 * std::max(std::max(sx,sy), sz) / 56) + 1; - const uint64_t psy = psx; + Vec3 normal(nx, ny, nz); + normal /= normal.norm(); Bbox2d bbx; - std::vector visited(psx * psy); + auto bases = create_orthonormal_basis(normal, anisotropy, positive_basis); + Vec3 basis1 = std::get<0>(bases); + Vec3 basis2 = std::get<1>(bases); - if (out == NULL) { - out = new LABEL[psx * psy](); - } + Bbox2d plane_bbx = compute_slice_plane(pos, basis1, basis2, sx, sy, sz); + + const int64_t psx = plane_bbx.sx(); + const int64_t psy = plane_bbx.sy(); + + LABEL* out = new LABEL[psx * psy](); if (px < 0 || px >= sx) { - return std::tuple(out, bbx); + return std::tuple(out, psx, bbx); } else if (py < 0 || py >= sy) { - return std::tuple(out, bbx); + return std::tuple(out, psx, bbx); } else if (pz < 0 || pz >= sz) { - return std::tuple(out, bbx); + return std::tuple(out, psx, bbx); } - const Vec3 pos(px, py, pz); - Vec3 normal(nx, ny, nz); - normal /= normal.norm(); - - auto bases = create_orthonormal_basis(normal, positive_basis); - Vec3 basis1 = std::get<0>(bases) * anisotropy; - Vec3 basis2 = std::get<1>(bases) * anisotropy; + bbx.x_min = psx - 1; + bbx.x_max = 0; + bbx.y_min = psy - 1; + bbx.y_max = 0; - uint64_t plane_pos_x = psx / 2; - uint64_t plane_pos_y = psy / 2; + for (int64_t y = 0; y < psy; y++) { + for (int64_t x = 0; x < psx; x++) { - bbx.x_min = plane_pos_x; - bbx.x_max = plane_pos_x; - bbx.y_min = plane_pos_y; - bbx.y_max = plane_pos_y; + float dx = static_cast(x) + static_cast(plane_bbx.x_min); + float dy = static_cast(y) + static_cast(plane_bbx.y_min); - uint64_t ploc = plane_pos_x + psx * plane_pos_y; + Vec3 cur = pos + basis1 * dx + basis2 * dy; - std::stack stack; - stack.push(ploc); - - 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 || cur.y < 0 || cur.z < 0) { - continue; - } - else if (cur.x >= sx || cur.y >= sy || cur.z >= sz) { - continue; - } - - bbx.x_min = std::min(bbx.x_min, static_cast(x)); - bbx.x_max = std::max(bbx.x_max, static_cast(x)); - bbx.y_min = std::min(bbx.y_min, static_cast(y)); - bbx.y_max = std::max(bbx.y_max, static_cast(y)); - - uint64_t loc = static_cast(cur.x) + sx * ( - static_cast(cur.y) + sy * static_cast(cur.z) - ); + if (cur.x < 0 || cur.y < 0 || cur.z < 0) { + continue; + } + else if (cur.x >= sx || cur.y >= sy || cur.z >= sz) { + continue; + } - out[ploc] = labels[loc]; + uint64_t loc = static_cast(cur.x) + sx * ( + static_cast(cur.y) + sy * static_cast(cur.z) + ); - uint64_t up = ploc - psx; - uint64_t down = ploc + psx; - uint64_t left = ploc - 1; - uint64_t right = ploc + 1; + bbx.x_min = std::min(bbx.x_min, x); + bbx.x_max = std::max(bbx.x_max, x); + bbx.y_min = std::min(bbx.y_min, y); + bbx.y_max = std::max(bbx.y_max, y); - 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); + uint64_t ploc = x + psx * y; + out[ploc] = labels[loc]; } } - return std::tuple(out, bbx); + return std::tuple(out, psx, bbx); } };