Skip to content

feat: add support for grids to cc3d #46

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
wants to merge 6 commits into
base: main
Choose a base branch
from
Open

feat: add support for grids to cc3d #46

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
aaf5cf8
feat: add support for grid size to connected_components
william-silversmith Jun 17, 2025
c4f65b4
fix: correct grid coloring for cc3d.connected_components
william-silversmith Jun 25, 2025
f4d3aba
test: check grid coloring works for cc3d
william-silversmith Jun 25, 2025
ca807f9
chore: whitespace
william-silversmith Jun 25, 2025
a7cdaa4
fix: correct estimate
william-silversmith Jun 25, 2025
2590100
feat: add grid support to color_connectivity_graph
william-silversmith Jun 25, 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
27 changes: 27 additions & 0 deletions automated_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -792,7 +792,34 @@ def test_spurious_branch_elimination():

assert np.all(recovered == arr)

def test_connected_component_grid():
import fastcrackle
labels = np.ones([20,20,3], dtype=np.uint32, order="F")

cc_labels, cc_per_slice, N = fastcrackle.connected_components(labels, 20, 20);
assert N == 3
assert cc_per_slice == [1,1,1]

cc_labels, cc_per_slice, N = fastcrackle.connected_components(labels, 10, 10);
assert N == 3 * 4
assert cc_per_slice == [4,4,4]

cc_labels, cc_per_slice, N = fastcrackle.connected_components(labels, 15, 15);
assert N == 3 * 4
assert cc_per_slice == [4,4,4]

cc_labels, cc_per_slice, N = fastcrackle.connected_components(labels, 1, 1);
assert N == 20 * 20 * 3

labels = np.ones([20,10,3], dtype=np.uint32, order="F")

cc_labels, cc_per_slice, N = fastcrackle.connected_components(labels, 10, 10);
assert N == 3 * 2
assert cc_per_slice == [2,2,2]

cc_labels, cc_per_slice, N = fastcrackle.connected_components(labels, 40, 40);
assert N == 3
assert cc_per_slice == [1,1,1]



Expand Down
180 changes: 121 additions & 59 deletions src/cc3d.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
* Author: William Silversmith
* Affiliation: Seung Lab, Princeton University
* Date: August 2018 - June 2019, 2021, Dec. 2022
* Dec. 2024, June 2025
*
* ----
* LICENSE
Expand Down Expand Up @@ -145,18 +146,26 @@ OUT* relabel(
}

uint64_t estimate_provisional_label_count_vcg(
const std::vector<uint8_t>& vcg, const int64_t sx
const std::vector<uint8_t>& vcg,
const int64_t sx, const int64_t sy, const int64_t sz,
const int64_t gx
) {
uint64_t count = 0; // number of transitions between labels
int64_t voxels = static_cast<int64_t>(vcg.size());

for (int64_t loc = 0; loc < voxels; loc += sx) {
for (int64_t loc = 0; loc < voxels; loc += gx) {
count += 1;
for (int64_t x = 1; x < sx; x++) {
for (int64_t x = 1; x < gx; x++) {
count += ((vcg[loc+x] & 0b0010) == 0);
}
}

// if gx overhangs sx, then need to add some corrections
if (sx % gx != 0) {
count += sy * sz;
count = std::min(count, static_cast<uint64_t>(voxels));
}

return count;
}

Expand All @@ -165,13 +174,23 @@ OUT* color_connectivity_graph(
const std::vector<uint8_t> &vcg, // voxel connectivity graph
const int64_t sx, const int64_t sy, const int64_t sz,
OUT* out_labels = NULL,
uint64_t &N = _dummy_N
uint64_t &N = _dummy_N,
int64_t gx = -1, int64_t gy = -1
) {

const int64_t sxy = sx * sy;
const int64_t voxels = sx * sy * sz;

uint64_t max_labels = estimate_provisional_label_count_vcg(vcg, sx) + 1;
gx = gx > 0 ? gx : sx;
gy = gy > 0 ? gy : sy;

gx = std::min(gx, sx);
gy = std::min(gy, sy);

const int64_t nx = (sx + gx - 1) / gx;
const int64_t ny = (sy + gy - 1) / gy;

uint64_t max_labels = 1 + estimate_provisional_label_count_vcg(vcg, sx, sy, sz, gx);
max_labels = std::min(max_labels, static_cast<uint64_t>(std::numeric_limits<OUT>::max()));

if (out_labels == NULL) {
Expand All @@ -180,39 +199,53 @@ OUT* color_connectivity_graph(

DisjointSet<OUT> equivalences(max_labels);

const int64_t B = -1;
const int64_t C = -sx;

OUT new_label = 0;
for (int64_t z = 0; z < sz; z++) {
new_label++;
equivalences.add(new_label);

for (int64_t x = 0; x < sx; x++) {
if (x > 0 && (vcg[x + sxy * z] & 0b0010) == 0) {
new_label++;
equivalences.add(new_label);
}
out_labels[x + sxy * z] = new_label;
}
for (int64_t gy_i = 0, sy_left = sy; gy_i < ny; gy_i++, sy_left -= gy) {
int64_t tail_y = std::min(gy, sy_left);

const int64_t B = -1;
const int64_t C = -sx;
for (int64_t gx_i = 0, sx_left = sx; gx_i < nx; gx_i++, sx_left -= gx) {
int64_t tail_x = std::min(gx, sx_left);

for (int64_t y = 1; y < sy; y++) {
for (int64_t x = 0; x < sx; x++) {
int64_t loc = x + sx * y + sxy * z;
new_label++;
equivalences.add(new_label);

if (x > 0 && (vcg[loc] & 0b0010)) {
out_labels[loc] = out_labels[loc+B];
if (y > 0 && (vcg[loc + C] & 0b0010) == 0 && (vcg[loc] & 0b1000)) {
equivalences.unify(out_labels[loc], out_labels[loc+C]);
for (int64_t x = 0; x < tail_x; x++) {
int64_t loc = x + gx_i * gx + sx * gy_i * gx + sxy * z;
if (x > 0 && (vcg[loc] & 0b0010) == 0) {
new_label++;
equivalences.add(new_label);
}
}
else if (y > 0 && vcg[loc] & 0b1000) {
out_labels[loc] = out_labels[loc+C];
}
else {
new_label++;
out_labels[loc] = new_label;
equivalences.add(new_label);
}
}

for (int64_t y = 1; y < tail_y; y++) {
for (int64_t gx_i = 0, sx_left = sx; gx_i < nx; gx_i++, sx_left -= gx) {
int64_t tail_x = std::min(gx, sx_left);

for (int64_t x = 0; x < tail_x; x++) {
int64_t loc = x + gx_i * gx + sx * (y + gy_i * gy) + sxy * z;

if (x > 0 && (vcg[loc] & 0b0010)) {
out_labels[loc] = out_labels[loc+B];
if (y > 0 && (vcg[loc + C] & 0b0010) == 0 && (vcg[loc] & 0b1000)) {
equivalences.unify(out_labels[loc], out_labels[loc+C]);
}
}
else if (y > 0 && vcg[loc] & 0b1000) {
out_labels[loc] = out_labels[loc+C];
}
else {
new_label++;
out_labels[loc] = new_label;
equivalences.add(new_label);
}
}
}
}
}
Expand All @@ -224,40 +257,59 @@ OUT* color_connectivity_graph(

template <typename LABEL>
uint64_t estimate_provisional_label_count(
const LABEL* in_labels, const int64_t sx, const int64_t voxels
const LABEL* in_labels,
const int64_t sx, const int64_t sy, const int64_t sz,
const int64_t gx
) {
uint64_t count = 0; // number of transitions between labels
const int64_t voxels = sx * sy * sz;

for (int64_t loc = 0; loc < voxels; loc += sx) {
for (int64_t loc = 0; loc < voxels; loc += gx) {
count += 1;
for (int64_t x = 1; x < sx; x++) {
for (int64_t x = 1; x < gx; x++) {
count += (in_labels[loc + x] != in_labels[loc + x - 1]);
}
}

// if gx overhangs sx, then need to add some corrections
if (sx % gx != 0) {
count += sy * sz;
count = std::min(count, static_cast<uint64_t>(voxels));
}

return count;
}

template <typename LABEL, typename OUT>
OUT* connected_components2d_4(
const LABEL* in_labels,
const int64_t sx, const int64_t sy, const int64_t sz,
OUT *out_labels = NULL,
const uint64_t start_label = 0, uint64_t &N = _dummy_N
) {
const LABEL* in_labels,
const int64_t sx, const int64_t sy, const int64_t sz,
OUT *out_labels = NULL,
const uint64_t start_label = 0, uint64_t &N = _dummy_N,
int64_t gx = -1, int64_t gy = -1
) {

const int64_t sxy = sx * sy;
const int64_t voxels = sx * sy * sz;

uint64_t max_labels = estimate_provisional_label_count<LABEL>(in_labels, sx, voxels) + 1;
gx = gx > 0 ? gx : sx;
gy = gy > 0 ? gy : sy;

gx = std::min(gx, sx);
gy = std::min(gy, sy);

const int64_t nx = (sx + gx - 1) / gx;
const int64_t ny = (sy + gy - 1) / gy;

uint64_t max_labels = 1 + estimate_provisional_label_count<LABEL>(in_labels, sx, sy, sz, gx);
max_labels = std::min(max_labels, static_cast<uint64_t>(std::numeric_limits<OUT>::max()));

DisjointSet<OUT> equivalences(max_labels);

if (out_labels == NULL) {
out_labels = new OUT[voxels]();
}

/*
Layout of forward pass mask.
A is the current location.
Expand All @@ -275,25 +327,33 @@ OUT* connected_components2d_4(

LABEL cur = 0;
for (int64_t z = 0; z < sz; z++) {
for (int64_t y = 0; y < sy; y++) {
for (int64_t x = 0; x < sx; x++) {
loc = x + sx * y + sxy * z;
cur = in_labels[loc];

if (x > 0 && cur == in_labels[loc + B]) {
out_labels[loc] = out_labels[loc + B];
if (y > 0 && cur == in_labels[loc + C] && cur != in_labels[loc + D]) {
equivalences.unify(out_labels[loc], out_labels[loc + C]);
for (int64_t gy_i = 0, sy_left = sy; gy_i < ny; gy_i++, sy_left -= gy) {
int64_t tail_y = std::min(gy, sy_left);

for (int64_t y = 0; y < tail_y; y++) {
for (int64_t gx_i = 0, sx_left = sx; gx_i < nx; gx_i++, sx_left -= gx) {
int64_t tail_x = std::min(sx_left, gx);

for (int64_t x = 0; x < tail_x; x++) {
loc = x + gx * gx_i + sx * (y + gy_i * gy) + sxy * z;
cur = in_labels[loc];

if (x > 0 && cur == in_labels[loc + B]) {
out_labels[loc] = out_labels[loc + B];
if (y > 0 && cur == in_labels[loc + C] && cur != in_labels[loc + D]) {
equivalences.unify(out_labels[loc], out_labels[loc + C]);
}
}
else if (y > 0 && cur == in_labels[loc + C]) {
out_labels[loc] = out_labels[loc + C];
}
else {
next_label++;
out_labels[loc] = next_label;
equivalences.add(out_labels[loc]);
}
}
}
else if (y > 0 && cur == in_labels[loc + C]) {
out_labels[loc] = out_labels[loc + C];
}
else {
next_label++;
out_labels[loc] = next_label;
equivalences.add(out_labels[loc]);
}
}
}
}
Expand All @@ -307,7 +367,8 @@ OUT* connected_components(
const int64_t sx, const int64_t sy, const int64_t sz,
std::vector<uint64_t> &num_components_per_slice,
OUT* out_labels = NULL,
uint64_t &N = _dummy_N
uint64_t &N = _dummy_N,
const int64_t gx = -1, const int64_t gy = -1
) {

const int64_t sxy = sx * sy;
Expand All @@ -323,7 +384,8 @@ OUT* connected_components(
(in_labels + sxy * z),
sx, sy, 1,
(out_labels + sxy * z),
N, tmp_N
N, tmp_N,
gx, gy
);
num_components_per_slice[z] = tmp_N;
N += tmp_N;
Expand Down
19 changes: 11 additions & 8 deletions src/crackle.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,8 @@ std::vector<unsigned char> compress_helper(
const bool optimize_pins = false,
const bool auto_bgcolor = true,
const bool manual_bgcolor = 0,
size_t parallel = 1
size_t parallel = 1,
int64_t grid_size = 2147483648
) {
const int64_t voxels = sx * sy * sz;

Expand All @@ -68,6 +69,8 @@ std::vector<unsigned char> compress_helper(
}
parallel = std::min(parallel, static_cast<size_t>(sz));

grid_size = pow(2, std::floor(log2(grid_size)));

CrackleHeader header(
/*format_version=*/CrackleHeader::current_version,
/*label_format=*/label_format,
Expand All @@ -79,12 +82,8 @@ std::vector<unsigned char> compress_helper(
/*sy=*/sy,
/*sz=*/sz,

// grid size is not yet supported, but will be
// used for within-slice random access.
// very large value so that way when we start
// using random access decoders, old formats
// will be backwards compatible (each slice is 1 grid).
/*grid_size*/2147483648,
// grid size must be a power of 2
/*grid_size*/grid_size,

/*num_label_bytes=*/0,
/*fortran_order*/fortran_order,
Expand Down Expand Up @@ -163,7 +162,11 @@ std::vector<unsigned char> compress_helper(
crack_crcs = std::move(crack_crcs_tmp);
}
else {
auto [labels_binary_tmp, crack_crcs_tmp] = crackle::labels::encode_flat<LABEL, STORED_LABEL>(labels, sx, sy, sz, parallel);
auto [labels_binary_tmp, crack_crcs_tmp] = crackle::labels::encode_flat<LABEL, STORED_LABEL>(
labels, sx, sy, sz,
header.grid_size,
parallel
);
labels_binary = std::move(labels_binary_tmp);
crack_crcs = std::move(crack_crcs_tmp);
}
Expand Down
Loading