Skip to content

perf(slice): reduce memory allocation + replace stack #7

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

Open
william-silversmith wants to merge 7 commits into
base: main
Choose a base branch
from
Open

perf(slice): reduce memory allocation + replace stack #7

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
e198554
wip: try looping instead of stack
william-silversmith Feb 23, 2024
22edec6
wip: experimenting with much smaller plane
william-silversmith Feb 24, 2024
dfd8456
wip: refactoring to transmit psx
william-silversmith Feb 24, 2024
00d5bc7
fix: warnings about uint64 comparisons
william-silversmith Feb 24, 2024
e132781
fix: remove print statements, add a margin for dx,dy
william-silversmith Feb 24, 2024
7979a13
fix: ensure we account for the diagnonal
william-silversmith Feb 24, 2024
3012a72
refactor: remove unused visited
william-silversmith Feb 24, 2024
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
32 changes: 7 additions & 25 deletions src/fastxs3d.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<uint64_t>(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<decltype(dtype), py::array::f_style>({ psx, psx });
auto out = reinterpret_cast<decltype(dtype)*>(arr.request().ptr);
auto data = reinterpret_cast<decltype(dtype)*>(labels.request().ptr);
std::fill(out, out + pvoxels, 0);

std::tuple<decltype(dtype)*, xs3d::Bbox2d> tup = xs3d::cross_section_projection<decltype(dtype)>(
auto tup = xs3d::cross_section_projection<decltype(dtype)>(
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++;

Expand All @@ -133,6 +113,8 @@ auto projection(
}
}

delete[] out;

return cutout.view(py::str(labels.dtype()));
};

Expand Down
233 changes: 119 additions & 114 deletions src/xs3d.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,10 +10,8 @@
#include <vector>
#include <stdexcept>

namespace {

static uint8_t _dummy_contact = false;

namespace xs3d {
class Vec3 {
public:
float x, y, z;
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -714,8 +720,13 @@ float* cross_section(
}

std::tuple<Vec3, Vec3> 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);
Expand Down Expand Up @@ -755,128 +766,122 @@ std::tuple<Vec3, Vec3> 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 <typename LABEL>
std::tuple<LABEL*, Bbox2d> cross_section_projection(
std::tuple<LABEL*, int64_t, Bbox2d> 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<uint64_t>(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<bool> 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<float>(x) + static_cast<float>(plane_bbx.x_min);
float dy = static_cast<float>(y) + static_cast<float>(plane_bbx.y_min);

uint64_t ploc = plane_pos_x + psx * plane_pos_y;
Vec3 cur = pos + basis1 * dx + basis2 * dy;

std::stack<uint64_t> 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<float>(x) - static_cast<float>(plane_pos_x);
float dy = static_cast<float>(y) - static_cast<float>(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<int64_t>(x));
bbx.x_max = std::max(bbx.x_max, static_cast<int64_t>(x));
bbx.y_min = std::min(bbx.y_min, static_cast<int64_t>(y));
bbx.y_max = std::max(bbx.y_max, static_cast<int64_t>(y));

uint64_t loc = static_cast<uint64_t>(cur.x) + sx * (
static_cast<uint64_t>(cur.y) + sy * static_cast<uint64_t>(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<uint64_t>(cur.x) + sx * (
static_cast<uint64_t>(cur.y) + sy * static_cast<uint64_t>(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);
}

};
Expand Down