From e19855490b84953d07c1e6584acf422c729f2ea9 Mon Sep 17 00:00:00 2001 From: William Silversmith Date: Fri, 23 Feb 2024 01:44:42 -0500 Subject: [PATCH 1/7] wip: try looping instead of stack --- src/xs3d.hpp | 75 +++++++++++++++------------------------------------- 1 file changed, 22 insertions(+), 53 deletions(-) diff --git a/src/xs3d.hpp b/src/xs3d.hpp index 70f9a05..d3d0cde 100644 --- a/src/xs3d.hpp +++ b/src/xs3d.hpp @@ -816,63 +816,32 @@ std::tuple cross_section_projection( bbx.y_min = plane_pos_y; bbx.y_max = plane_pos_y; - uint64_t ploc = plane_pos_x + psx * plane_pos_y; - - 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; - } + for (float y = 0; y < psy; y++) { + for (float x = 0; x < psx; x++) { - 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)); + float dx = static_cast(x) - static_cast(plane_pos_x); + float dy = static_cast(y) - static_cast(plane_pos_y); - uint64_t loc = static_cast(cur.x) + sx * ( - static_cast(cur.y) + sy * static_cast(cur.z) - ); + Vec3 cur = pos + basis1 * dx + basis2 * dy; - out[ploc] = labels[loc]; - - uint64_t up = ploc - psx; - uint64_t down = ploc + psx; - uint64_t left = ploc - 1; - uint64_t right = ploc + 1; + if (cur.x < 0 || cur.y < 0 || cur.z < 0) { + continue; + } + else if (cur.x >= sx || cur.y >= sy || cur.z >= sz) { + continue; + } - 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); + 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) + ); + uint64_t ploc = x + psx * y; + + out[ploc] = labels[loc]; } } From 22edec683d194dcd5f4a31397fdcf62c9557cf3b Mon Sep 17 00:00:00 2001 From: William Silversmith Date: Fri, 23 Feb 2024 20:41:12 -0500 Subject: [PATCH 2/7] wip: experimenting with much smaller plane --- src/fastxs3d.cpp | 55 ++++++++------- src/xs3d.hpp | 169 ++++++++++++++++++++++++++++++----------------- 2 files changed, 138 insertions(+), 86 deletions(-) diff --git a/src/fastxs3d.cpp b/src/fastxs3d.cpp index 59165f5..5029a9a 100644 --- a/src/fastxs3d.cpp +++ b/src/fastxs3d.cpp @@ -83,20 +83,25 @@ auto projection( ? 1 : labels.shape()[2]; - float wx = anisotropy.at(0); - float wy = anisotropy.at(1); - float wz = anisotropy.at(2); + xs3d::Vec3 aniso(anisotropy.at(0), anisotropy.at(1), anisotropy.at(2)); + xs3d::Vec3 pos(point.at(0), point.at(1), point.at(2)); + xs3d::Vec3 norm(normal.at(0), normal.at(1), normal.at(2)); + norm /= norm.norm(); + + std::tuple bases = xs3d::create_orthonormal_basis(norm, aniso, standardize_basis); + xs3d::Bbox2d plane_bbx = xs3d::compute_slice_plane( + pos, + std::get<0>(bases), std::get<1>(bases), + sx, sy, sz + ); - 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 int64_t psx = plane_bbx.sx(); + const int64_t psy = plane_bbx.sy(); + const int64_t pvoxels = psx * psy; - const uint64_t distortion = static_cast(ceil(maxval)); + plane_bbx.print(); - // 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; + printf("fastxs3d.cpp psx %d psy %d\n", psx, psy); py::array arr; @@ -116,24 +121,24 @@ auto projection( out ); - xs3d::Bbox2d bbox = std::get<1>(tup); - bbox.x_max++; - bbox.y_max++; + // xs3d::Bbox2d bbox = std::get<1>(tup); + // bbox.x_max++; + // bbox.y_max++; - auto cutout = py::array_t({ bbox.sx(), bbox.sy() }); - auto cutout_ptr = reinterpret_cast(cutout.request().ptr); + // auto cutout = py::array_t({ bbox.sx(), bbox.sy() }); + // auto cutout_ptr = reinterpret_cast(cutout.request().ptr); - int64_t csx = bbox.sx(); + // int64_t csx = bbox.sx(); - for (int64_t y = bbox.y_min; y < bbox.y_max; y++) { - for (int64_t x = bbox.x_min; x < bbox.x_max; x++) { - cutout_ptr[ - (x - bbox.x_min) + csx * (y - bbox.y_min) - ] = out[x + psx * y]; - } - } + // for (int64_t y = bbox.y_min; y < bbox.y_max; y++) { + // for (int64_t x = bbox.x_min; x < bbox.x_max; x++) { + // cutout_ptr[ + // (x - bbox.x_min) + csx * (y - bbox.y_min) + // ] = out[x + psx * y]; + // } + // } - return cutout.view(py::str(labels.dtype())); + return arr.view(py::str(labels.dtype())); }; int data_width = labels.dtype().itemsize(); diff --git a/src/xs3d.hpp b/src/xs3d.hpp index d3d0cde..a17bbd4 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,7 +766,46 @@ 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 / basis.x; + float dxy = -pos.y / basis.y; + float dxz = -pos.z / 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 = (sx-pos.x) / basis.x; + float dyy = (sy-pos.y) / basis.y; + float dyz = (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); + float dx_max = maxfn(basis1); + + float dy_min = minfn(basis2); + float dy_max = maxfn(basis2); + + return Bbox2d(dx_min, dx_max, dy_min, dy_max); } template @@ -769,22 +819,22 @@ std::tuple cross_section_projection( const bool positive_basis, LABEL* out = NULL ) { + 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); + + Bbox2d plane_bbx = compute_slice_plane(pos, basis1, basis2, sx, sy, sz); + + const uint64_t psx = plane_bbx.sx(); + const uint64_t psy = plane_bbx.sy(); if (out == NULL) { out = new LABEL[psx * psy](); @@ -800,30 +850,25 @@ std::tuple cross_section_projection( return std::tuple(out, 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; - - uint64_t plane_pos_x = psx / 2; - uint64_t plane_pos_y = psy / 2; + std::vector visited(psx * psy); - bbx.x_min = plane_pos_x; - bbx.x_max = plane_pos_x; - bbx.y_min = plane_pos_y; - bbx.y_max = plane_pos_y; + bbx.x_min = sx; + bbx.x_max = 0; + bbx.y_min = sy; + bbx.y_max = 0; - for (float y = 0; y < psy; y++) { - for (float x = 0; x < psx; x++) { + // plane_bbx.print(); + // pos.print("pos"); - float dx = static_cast(x) - static_cast(plane_pos_x); - float dy = static_cast(y) - static_cast(plane_pos_y); + for (int64_t y = 0; y < psy; y++) { + for (int64_t x = 0; x < psx; x++) { - Vec3 cur = pos + basis1 * dx + basis2 * dy; + float dx = static_cast(x) + static_cast(plane_bbx.x_min); + float dy = static_cast(y) + static_cast(plane_bbx.y_min); + Vec3 cur = pos + basis1 * dx + basis2 * dy; + // printf("%.1f %.1f\n", dx, dy); + // cur.print("cur"); if (cur.x < 0 || cur.y < 0 || cur.z < 0) { continue; } @@ -831,10 +876,10 @@ std::tuple cross_section_projection( 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)); + 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); uint64_t loc = static_cast(cur.x) + sx * ( static_cast(cur.y) + sy * static_cast(cur.z) @@ -845,6 +890,8 @@ std::tuple cross_section_projection( } } + // bbx.print(); + return std::tuple(out, bbx); } From dfd845614fe7db964469b2f3fbfa4b771a9ed0b0 Mon Sep 17 00:00:00 2001 From: William Silversmith Date: Fri, 23 Feb 2024 21:07:37 -0500 Subject: [PATCH 3/7] wip: refactoring to transmit psx --- src/fastxs3d.cpp | 63 +++++++++++++++--------------------------------- src/xs3d.hpp | 18 ++++++-------- 2 files changed, 28 insertions(+), 53 deletions(-) diff --git a/src/fastxs3d.cpp b/src/fastxs3d.cpp index 5029a9a..7e8416d 100644 --- a/src/fastxs3d.cpp +++ b/src/fastxs3d.cpp @@ -83,62 +83,39 @@ auto projection( ? 1 : labels.shape()[2]; - xs3d::Vec3 aniso(anisotropy.at(0), anisotropy.at(1), anisotropy.at(2)); - xs3d::Vec3 pos(point.at(0), point.at(1), point.at(2)); - xs3d::Vec3 norm(normal.at(0), normal.at(1), normal.at(2)); - norm /= norm.norm(); - - std::tuple bases = xs3d::create_orthonormal_basis(norm, aniso, standardize_basis); - xs3d::Bbox2d plane_bbx = xs3d::compute_slice_plane( - pos, - std::get<0>(bases), std::get<1>(bases), - sx, sy, sz - ); - - const int64_t psx = plane_bbx.sx(); - const int64_t psy = plane_bbx.sy(); - const int64_t pvoxels = psx * psy; - - plane_bbx.print(); - - printf("fastxs3d.cpp psx %d psy %d\n", psx, psy); - - 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); - // bbox.x_max++; - // bbox.y_max++; + 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++; - // auto cutout = py::array_t({ bbox.sx(), bbox.sy() }); - // auto cutout_ptr = reinterpret_cast(cutout.request().ptr); + auto cutout = py::array_t({ bbox.sx(), bbox.sy() }); + auto cutout_ptr = reinterpret_cast(cutout.request().ptr); - // int64_t csx = bbox.sx(); + int64_t csx = bbox.sx(); - // for (int64_t y = bbox.y_min; y < bbox.y_max; y++) { - // for (int64_t x = bbox.x_min; x < bbox.x_max; x++) { - // cutout_ptr[ - // (x - bbox.x_min) + csx * (y - bbox.y_min) - // ] = out[x + psx * y]; - // } - // } + for (int64_t y = bbox.y_min; y < bbox.y_max; y++) { + for (int64_t x = bbox.x_min; x < bbox.x_max; x++) { + cutout_ptr[ + (x - bbox.x_min) + csx * (y - bbox.y_min) + ] = out[x + psx * y]; + } + } - return arr.view(py::str(labels.dtype())); + delete[] out; + + return cutout.view(py::str(labels.dtype())); }; int data_width = labels.dtype().itemsize(); diff --git a/src/xs3d.hpp b/src/xs3d.hpp index a17bbd4..f59abfe 100644 --- a/src/xs3d.hpp +++ b/src/xs3d.hpp @@ -809,15 +809,14 @@ Bbox2d compute_slice_plane( } 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); @@ -836,18 +835,16 @@ std::tuple cross_section_projection( const uint64_t psx = plane_bbx.sx(); const uint64_t psy = plane_bbx.sy(); - if (out == NULL) { - out = new LABEL[psx * psy](); - } + 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); } std::vector visited(psx * psy); @@ -884,6 +881,7 @@ std::tuple cross_section_projection( uint64_t loc = static_cast(cur.x) + sx * ( static_cast(cur.y) + sy * static_cast(cur.z) ); + // printf("%d, %d\n", x,y); uint64_t ploc = x + psx * y; out[ploc] = labels[loc]; @@ -892,7 +890,7 @@ std::tuple cross_section_projection( // bbx.print(); - return std::tuple(out, bbx); + return std::tuple(out, psx, bbx); } }; From 00d5bc771d7504f40f88253c911a56b616005acb Mon Sep 17 00:00:00 2001 From: William Silversmith Date: Fri, 23 Feb 2024 21:11:17 -0500 Subject: [PATCH 4/7] fix: warnings about uint64 comparisons --- src/xs3d.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/xs3d.hpp b/src/xs3d.hpp index f59abfe..e9b8517 100644 --- a/src/xs3d.hpp +++ b/src/xs3d.hpp @@ -809,7 +809,7 @@ Bbox2d compute_slice_plane( } 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, @@ -832,8 +832,8 @@ std::tuple cross_section_projection( Bbox2d plane_bbx = compute_slice_plane(pos, basis1, basis2, sx, sy, sz); - const uint64_t psx = plane_bbx.sx(); - const uint64_t psy = plane_bbx.sy(); + const int64_t psx = plane_bbx.sx(); + const int64_t psy = plane_bbx.sy(); LABEL* out = new LABEL[psx * psy](); From e1327817b36a5a373ca88298afab9c0fd326c8ce Mon Sep 17 00:00:00 2001 From: William Silversmith Date: Fri, 23 Feb 2024 22:19:53 -0500 Subject: [PATCH 5/7] fix: remove print statements, add a margin for dx,dy --- src/xs3d.hpp | 29 ++++++++++++----------------- 1 file changed, 12 insertions(+), 17 deletions(-) diff --git a/src/xs3d.hpp b/src/xs3d.hpp index e9b8517..69b88ad 100644 --- a/src/xs3d.hpp +++ b/src/xs3d.hpp @@ -799,11 +799,11 @@ Bbox2d compute_slice_plane( return ceil(std::max(std::max(dyx, dyy), dyz)); }; - float dx_min = minfn(basis1); - float dx_max = maxfn(basis1); + float dx_min = minfn(basis1) - 1; + float dx_max = maxfn(basis1) + 1; - float dy_min = minfn(basis2); - float dy_max = maxfn(basis2); + float dy_min = minfn(basis2) - 1; + float dy_max = maxfn(basis2) + 1; return Bbox2d(dx_min, dx_max, dy_min, dy_max); } @@ -849,14 +849,11 @@ std::tuple cross_section_projection( std::vector visited(psx * psy); - bbx.x_min = sx; + bbx.x_min = psx - 1; bbx.x_max = 0; - bbx.y_min = sy; + bbx.y_min = psy - 1; bbx.y_max = 0; - // plane_bbx.print(); - // pos.print("pos"); - for (int64_t y = 0; y < psy; y++) { for (int64_t x = 0; x < psx; x++) { @@ -864,8 +861,7 @@ std::tuple cross_section_projection( float dy = static_cast(y) + static_cast(plane_bbx.y_min); Vec3 cur = pos + basis1 * dx + basis2 * dy; - // printf("%.1f %.1f\n", dx, dy); - // cur.print("cur"); + if (cur.x < 0 || cur.y < 0 || cur.z < 0) { continue; } @@ -873,22 +869,21 @@ std::tuple cross_section_projection( continue; } + uint64_t loc = static_cast(cur.x) + sx * ( + static_cast(cur.y) + sy * static_cast(cur.z) + ); + 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); - uint64_t loc = static_cast(cur.x) + sx * ( - static_cast(cur.y) + sy * static_cast(cur.z) - ); - // printf("%d, %d\n", x,y); uint64_t ploc = x + psx * y; - out[ploc] = labels[loc]; } } - // bbx.print(); + bbx.print(); return std::tuple(out, psx, bbx); } From 7979a139f03ec675fae2fe82c00bf4f64272665c Mon Sep 17 00:00:00 2001 From: William Silversmith Date: Fri, 23 Feb 2024 23:33:51 -0500 Subject: [PATCH 6/7] fix: ensure we account for the diagnonal --- src/xs3d.hpp | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/xs3d.hpp b/src/xs3d.hpp index 69b88ad..7887585 100644 --- a/src/xs3d.hpp +++ b/src/xs3d.hpp @@ -776,9 +776,9 @@ Bbox2d compute_slice_plane( ) { auto minfn = [&](const Vec3& basis) { - float dxx = -pos.x / basis.x; - float dxy = -pos.y / basis.y; - float dxz = -pos.z / basis.z; + 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; @@ -788,9 +788,9 @@ Bbox2d compute_slice_plane( }; auto maxfn = [&](const Vec3& basis) { - float dyx = (sx-pos.x) / basis.x; - float dyy = (sy-pos.y) / basis.y; - float dyz = (sz-pos.z) / basis.z; + 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; @@ -883,8 +883,6 @@ std::tuple cross_section_projection( } } - bbx.print(); - return std::tuple(out, psx, bbx); } From 3012a72f4b3a4e763c9f9ad9375d3fa76324d21b Mon Sep 17 00:00:00 2001 From: William Silversmith Date: Fri, 23 Feb 2024 23:39:05 -0500 Subject: [PATCH 7/7] refactor: remove unused visited --- src/xs3d.hpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/xs3d.hpp b/src/xs3d.hpp index 7887585..d78fcdc 100644 --- a/src/xs3d.hpp +++ b/src/xs3d.hpp @@ -847,8 +847,6 @@ std::tuple cross_section_projection( return std::tuple(out, psx, bbx); } - std::vector visited(psx * psy); - bbx.x_min = psx - 1; bbx.x_max = 0; bbx.y_min = psy - 1;