seung-lab / euclidean-distance-transform-3d

Euclidean distance & signed distance transform for multi-label 3D anisotropic images using marching parabolas.
GNU General Public License v3.0
232 stars 37 forks source link

Question: reproducing Matlab's bwdist #55

Closed DoubleEspresso closed 1 month ago

DoubleEspresso commented 1 month ago

Hello,

I've tried using your edtsq on uint8 input to reproduce the output of Matlab's bwdist. Bwdist will compute the min distance between background (label 0) and nearest mask pt (label 1). In otherwords the output of sqrt( edtsq(mask) ) = Matlab::bwdist(~mask).

My hope was I could take the complement of input mask and edtsq would reproduce this distance map, but that doesn't seem to work, adding a border also does not work unless I did it incorrectly.

Is there a simple edit to the code/method that could be added to support this?

Thanks in advance!

william-silversmith commented 1 month ago

Hi! After reading the Matlab documentation, I think this does the trick:

edt.edt(mask == 0)

Let me know if you need additional help.

DoubleEspresso commented 1 month ago

Thanks for the quick response!

I did actually try this - I've attached an image for reference. Image on the left is output from bwdist, while image on the right is the result of edt.edt(mask == 0). White pixels in the right image are all inf.

mask

william-silversmith commented 1 month ago

Can you send me the mask you're using? I tried it locally with the example on the Matlab website and it worked. I can probably tell you what's up if I can reproduce this locally.

DoubleEspresso commented 1 month ago

Sure - attached is the small test volume I'm using.

I've checked out 2.4.1 in CPP for this too, FYI.

mask.zip

william-silversmith commented 1 month ago

I tried this with mask18.tiff, and got this result (the purple area is the mask overlaid).

image

I got this after converting mask18.tiff to mask18.png using imagemagick. Here's the code:

import numpy as np
import edt
import pyspng
import microviewer

with open("mask/mask18.png", "rb") as f:
    mask = pyspng.load(f.read()).astype(bool).T

dt = edt.edt(mask == 0)
microviewer.hyperview(dt, mask)

To reproduce:

pip install numpy edt pyspng-seunglab microviewer

Is this what you wanted? Am I doing something different?

DoubleEspresso commented 1 month ago

This looks good! I must be messing something up - here's the snippet I'm using in C++, is this wrong?

dst.Init(volSize, 0);
auto inMask = const_cast<uint8*>(src.ptr());
auto output = const_cast<float*>(dst.ptr());

edt::edtsq<uint8>(inMask, dst.width, dst.height, dst.depth, 1.0, 1.0, 1.0, false, 1, output);

dst.Sqrt(); 
william-silversmith commented 1 month ago

Ah doing it in C++ adds another few wrinkles. I should probably note that the cpp code is basically never tested, but the python C++ version is. I should probably delete the "C++" version as it is likely to mislead people. Grab python/edt.hpp as it is most likely error free and contains all the same functionality.

One thing to keep in mind is that the function expects the data in Fortran order. If your data is in C order you need to reverse the dimensions and anisotropy (by coincidence, the algorithm is fine with that). Since your data is also binary, you can benefit from a small perf increase via binary_edt (I probably should just make that a template specialization for bool* — this was one of my first C++ projects).

DoubleEspresso commented 1 month ago

Excellent, thank you! I'll checkout python/edt.hpp and let you know how it goes.

DoubleEspresso commented 1 month ago

Looks perfect! Checked both the edt and edtsq methods. Might want to add a note in the readme for people looking for the CPP codes directly. Thanks again!

william-silversmith commented 1 month ago

Fabulous! I'll have to do some house cleaning and edit the README.

william-silversmith commented 1 month ago

done!