Skip to content
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

Improve performance of binary EDT on black pixels #6

Merged
merged 3 commits into from
Sep 11, 2018
Merged

Improve performance of binary EDT on black pixels #6

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
85e4b4b
perf: efficiently deal with runs of black on binary images
william-silversmith Sep 11, 2018
09d167b
perf: accelerate binary EDT on starting black pixels
william-silversmith Sep 11, 2018
c4cc51d
test: add binary EDT to scipy comparisons
william-silversmith Sep 11, 2018
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
63 changes: 43 additions & 20 deletions cpp/edt.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -460,9 +460,11 @@ float* _binary_edt3dsq(T* binaryimg,
const size_t sxy = sx * sy;
const size_t voxels = sz * sxy;

size_t x,y,z;

float *workspace = new float[sx * sy * sz]();
for (size_t z = 0; z < sz; z++) {
for (size_t y = 0; y < sy; y++) {
for (z = 0; z < sz; z++) {
for (y = 0; y < sy; y++) {
// Might be possible to write this as a single pass, might be faster
// however, it's already only using about 3-5% of total CPU time.
// NOTE: Tried it, same speed overall.
Expand All @@ -478,24 +480,38 @@ float* _binary_edt3dsq(T* binaryimg,
tofinite(workspace, voxels);
}

for (size_t z = 0; z < sz; z++) {
for (size_t x = 0; x < sx; x++) {
size_t offset;
for (z = 0; z < sz; z++) {
for (x = 0; x < sx; x++) {
offset = x + sxy * z;
for (y = 0; y < sy; y++) {
if (workspace[offset + sx*y]) {
break;
}
}

_squared_edt_1d_parabolic(
(workspace + x + sxy * z),
(workspace + x + sxy * z),
sy, sx, wy,
black_border, black_border
(workspace + offset + sx * y),
(workspace + offset + sx * y),
sy - y, sx, wy,
black_border || (y > 0), black_border
);
}
}

for (size_t y = 0; y < sy; y++) {
for (size_t x = 0; x < sx; x++) {
for (y = 0; y < sy; y++) {
for (x = 0; x < sx; x++) {
offset = x + sx * y;
for (z = 0; z < sz; z++) {
if (workspace[offset + sxy*z]) {
break;
}
}
_squared_edt_1d_parabolic(
(workspace + x + sx * y),
(workspace + x + sx * y),
sz, sxy, wz,
black_border, black_border
(workspace + offset + sxy * z),
(workspace + offset + sxy * z),
sz - z, sxy, wz,
black_border || (z > 0), black_border
);
}
}
Expand Down Expand Up @@ -598,9 +614,10 @@ float* _binary_edt2dsq(T* binaryimg,
const bool black_border=false) {

const size_t voxels = sx * sy;
size_t x,y;

float *workspace = new float[sx * sy]();
for (size_t y = 0; y < sy; y++) {
for (y = 0; y < sy; y++) {
squared_edt_1d_multi_seg<T>(
(binaryimg + sx * y), (workspace + sx * y),
sx, 1, wx, black_border
Expand All @@ -611,12 +628,18 @@ float* _binary_edt2dsq(T* binaryimg,
tofinite(workspace, voxels);
}

for (size_t x = 0; x < sx; x++) {
for (x = 0; x < sx; x++) {
for (y = 0; y < sy; y++) {
if (workspace[x + y * sx]) {
break;
}
}

_squared_edt_1d_parabolic(
(workspace + x),
(workspace + x),
sy, sx, wy,
black_border, black_border
(workspace + x + y * sx),
(workspace + x + y * sx),
sy - y, sx, wy,
black_border || (y > 0), black_border
);
}

Expand Down
63 changes: 43 additions & 20 deletions python/edt.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -460,9 +460,11 @@ float* _binary_edt3dsq(T* binaryimg,
const size_t sxy = sx * sy;
const size_t voxels = sz * sxy;

size_t x,y,z;

float *workspace = new float[sx * sy * sz]();
for (size_t z = 0; z < sz; z++) {
for (size_t y = 0; y < sy; y++) {
for (z = 0; z < sz; z++) {
for (y = 0; y < sy; y++) {
// Might be possible to write this as a single pass, might be faster
// however, it's already only using about 3-5% of total CPU time.
// NOTE: Tried it, same speed overall.
Expand All @@ -478,24 +480,38 @@ float* _binary_edt3dsq(T* binaryimg,
tofinite(workspace, voxels);
}

for (size_t z = 0; z < sz; z++) {
for (size_t x = 0; x < sx; x++) {
size_t offset;
for (z = 0; z < sz; z++) {
for (x = 0; x < sx; x++) {
offset = x + sxy * z;
for (y = 0; y < sy; y++) {
if (workspace[offset + sx*y]) {
break;
}
}

_squared_edt_1d_parabolic(
(workspace + x + sxy * z),
(workspace + x + sxy * z),
sy, sx, wy,
black_border, black_border
(workspace + offset + sx * y),
(workspace + offset + sx * y),
sy - y, sx, wy,
black_border || (y > 0), black_border
);
}
}

for (size_t y = 0; y < sy; y++) {
for (size_t x = 0; x < sx; x++) {
for (y = 0; y < sy; y++) {
for (x = 0; x < sx; x++) {
offset = x + sx * y;
for (z = 0; z < sz; z++) {
if (workspace[offset + sxy*z]) {
break;
}
}
_squared_edt_1d_parabolic(
(workspace + x + sx * y),
(workspace + x + sx * y),
sz, sxy, wz,
black_border, black_border
(workspace + offset + sxy * z),
(workspace + offset + sxy * z),
sz - z, sxy, wz,
black_border || (z > 0), black_border
);
}
}
Expand Down Expand Up @@ -598,9 +614,10 @@ float* _binary_edt2dsq(T* binaryimg,
const bool black_border=false) {

const size_t voxels = sx * sy;
size_t x,y;

float *workspace = new float[sx * sy]();
for (size_t y = 0; y < sy; y++) {
for (y = 0; y < sy; y++) {
squared_edt_1d_multi_seg<T>(
(binaryimg + sx * y), (workspace + sx * y),
sx, 1, wx, black_border
Expand All @@ -611,12 +628,18 @@ float* _binary_edt2dsq(T* binaryimg,
tofinite(workspace, voxels);
}

for (size_t x = 0; x < sx; x++) {
for (x = 0; x < sx; x++) {
for (y = 0; y < sy; y++) {
if (workspace[x + y * sx]) {
break;
}
}

_squared_edt_1d_parabolic(
(workspace + x),
(workspace + x),
sy, sx, wy,
black_border, black_border
(workspace + x + y * sx),
(workspace + x + y * sx),
sy - y, sx, wy,
black_border || (y > 0), black_border
);
}

Expand Down
10 changes: 6 additions & 4 deletions python/perf.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,15 +5,17 @@

import edt

labels = np.ones(shape=(512, 512, 512), dtype=np.uint32)
labels = np.ones(shape=(4096, 4096), dtype=np.uint32)

start = time.time()
res = edt.edtsq(labels, anisotropy=(1,1,1))
print('Multi-label EDT: ', time.time() - start, ' sec.')

for i in range(300):
x = (labels == i) * res
binlabels = labels.astype(np.bool)

print('Multi-label EDT: ', time.time() - start, ' sec.')
start = time.time()
res = edt.edtsq(binlabels, anisotropy=(1,1,1))
print('Binary EDT: ', time.time() - start, ' sec.')

start = time.time()
ndimage.distance_transform_edt(labels)
Expand Down
58 changes: 31 additions & 27 deletions python/test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -399,21 +399,22 @@ def test_2d_scipy_comparison_black_border():

def test_2d_scipy_comparison():
for _ in range(20):
randos = np.random.randint(0, 2, size=(5, 5), dtype=np.uint32)
labels = np.zeros( (randos.shape[0] + 2, randos.shape[1] + 2), dtype=np.uint32)
for dtype in (np.uint32, np.bool):
randos = np.random.randint(0, 2, size=(5, 5), dtype=dtype)
labels = np.zeros( (randos.shape[0] + 2, randos.shape[1] + 2), dtype=dtype)

print("INPUT")
print(labels)
print("INPUT")
print(labels)

print("MLAEDT")
mlaedt_result = edt.edt(labels, black_border=False)
print(mlaedt_result)
print("MLAEDT")
mlaedt_result = edt.edt(labels, black_border=False)
print(mlaedt_result)

print("SCIPY")
scipy_result = ndimage.distance_transform_edt(labels)
print(scipy_result)
print("SCIPY")
scipy_result = ndimage.distance_transform_edt(labels)
print(scipy_result)

assert np.all( np.abs(scipy_result - mlaedt_result) < 0.000001 )
assert np.all( np.abs(scipy_result - mlaedt_result) < 0.000001 )

def test_three_d():
def cmp(labels, ans, types=TYPES, anisotropy=(1.0, 1.0, 1.0)):
Expand Down Expand Up @@ -536,24 +537,27 @@ def cmp(labels, ans, types=TYPES, anisotropy=(1.0, 1.0, 1.0)):
], anisotropy=(4,40,4))

def test_3d_scipy_comparison():
randos = np.random.randint(0, 2, size=(10, 10, 10), dtype=np.uint32)
labels = np.zeros( (randos.shape[0] + 2, randos.shape[1] + 2, randos.shape[2] + 2), dtype=np.uint32)
# Scipy requires zero borders
labels[1:-1,1:-1,1:-1] = randos
for _ in range(20):
for dtype in (np.uint32, np.bool):
randos = np.random.randint(0, 2, size=(10, 10, 10), dtype=dtype)
labels = np.zeros( (randos.shape[0] + 2, randos.shape[1] + 2, randos.shape[2] + 2), dtype=dtype)
# Scipy requires zero borders
labels[1:-1,1:-1,1:-1] = randos

print("INPUT")
print(labels)

print("INPUT")
print(labels)
print("MLAEDT")
mlaedt_result = edt.edt(labels, black_border=False)
print(mlaedt_result)

print("MLAEDT")
mlaedt_result = edt.edt(labels, black_border=False)
print(mlaedt_result)
print("SCIPY")
scipy_result = ndimage.distance_transform_edt(labels)
print(scipy_result)

print("SCIPY")
scipy_result = ndimage.distance_transform_edt(labels)
print(scipy_result)
print("DIFF")
print(np.abs(scipy_result == mlaedt_result))
print(np.max(np.abs(scipy_result - mlaedt_result)))

print("DIFF")
print(np.abs(scipy_result == mlaedt_result))
print(np.max(np.abs(scipy_result - mlaedt_result)))
assert np.all( np.abs(scipy_result - mlaedt_result) < 0.000001 )

assert np.all( np.abs(scipy_result - mlaedt_result) < 0.000001 )