`utils/dt_morphometry` refactoring? Delete?

[njneeteson]

are these functions used elsewhere in the project? I have a new (much faster) function for calculating local thickness and most of the other functions in that file should probably be in bone_morphometry and not in the utils module.

[mkuczyns]

The calc_mean_thickness function is used in the joint_space_analysis module which is why I added those functions to the utils module. I've pushed the updated package structure to GitHub so feel free to upload your new functions or move things around now.

If we decide to move the bone morphometry functions out of the utils module, we'll need to update the imports in the jsw_morphometry.py script and Jupyter Notebook in examples/jsw_notebook.ipynb.

[njneeteson]

OK yeah since it is used all over probably makes more sense to stay in utils.

Do you mind if I completely replace everything that's in there with my new code though? I have something that's a lot faster (necessary for morphometry to be feasible on radius/tibia or anything larger) and that works with non-isotropic resolution images. I'll make sure to refactor so that anything that uses it still works.

[njneeteson]

Actually now that I'm looking through the code you already have a duplicate dt_morphometry sub-module in the joint_space_analysis module that those joint space width calculating sub-modules import the mean thickness calculation functions from (rather than from utils/...). The code in utils/dt_morphometry is not used anywhere in the package other than being imported in an __init__ somewhere. Do you want to leave this as-is or should the same distance transform logic be used everywhere?

[mkuczyns]

That's my mistake! There shouldn't be another copy of dt_morphometry in the joint_space_analysis module. You can upload your changes to util. I've pushed a new commit to remove the extra copy of dt_morphometry and updated the imports for the JSW module.

[njneeteson]

No problem.. OK all of my changes are good to be uploaded once I have access to the repository so I can push my branch to the origin here. Then you can take a look and see if you want to merge. There's a lot of changes and you might want to test your stuff to make sure it's not too affected before merging.

[mkuczyns]

Just added you as an admin!

[njneeteson]

Thanks! I pushed up my branch. Take a look when you have time, no rush

[njneeteson]

By the way I was looking at your notebook to see if what I had changed had changed you results. Looks like they did.. the mean thickness values reported were correct for a range of cylinders and spheres so I'm not sure what's up but it needs some debugging / investigation before it can be merged.

[njneeteson]

OK I changed a bunch of stuff and it's better now. The notebook / JSW analysis still is messed up and I'm not sure why, but I did manage to generate this nice image applying it to one of our lab test aims:

So I'm not sure what's up with the notebook and why my changes have messed it up.. I'll keep looking at it tomorrow.

[njneeteson]

OK I fiddled a lot more and it's again much better. Tested against synthetic spheres, cylinders, and plates it pretty much always gets the average radius/thickness correct. I also loaded up your JC_MASK.mha your notebook creates and tested it on that:

Still the notebook seems like I've messed it up. I might need you to take a look and see if you can figure out what's wrong. Or maybe use my branch to run your JSW analysis on an image or two and see if it's messed up your stats. At this point I can't see what's wrong with the local thickness calculation so I've got no idea what the problem is. I tried to simplify the calculation of the stats so it's possible that the way the mean/std/min/max has been messed up and maybe you can correct that so it's consistent with how it was before.

[mkuczyns]

Looks like of the JSW parameters have changed with your branch. Here is a comparison of the IPL method, the "old" Python method, and "new" method using your branch. I'll take a look at your branch to see what's changed.

[njneeteson]

Yeah.. I am pretty sure the local thickness / Hildebrand algorithm part is all good and I probably just messed up how the mean/std/min/max are calculated afterwards. It was so convoluted in the old Python implementation and I tried to simplify it but I guess I made some mistake somewhere. Maybe you can spot it

[njneeteson]

I haven't tested the code with a full microarchitectural analysis but I did just calculate Ct.Th and Tb.Th for an HR-pQCT image (from NORMXTII), with the new code:

           |  Python  | Reference (IPL)
----------------------------------------
Ct.Th [mm] |    1.395 | 1.379 
Tb.Th [mm] |    0.281 | 0.271 

That is really close and kind of leads me to believe the new code is outputting at least the correct mean thickness for the mask it is given. Is it possible there was an error in the old version and now it is more accurate? Though that doesn't explain why it was consistent with IPL JSW calcs and isn't anymore.

[njneeteson]

One of the bigger things I changed from the old version to the new version is the skeletonization algorithm. It used to be done using skeletonize_3d from skimage and I swapped that out for BinaryThinning from SimpleITK. The main reason I made this swap is because the skeletonizations for flat plates with the skimage version looked really poor and BinaryThinning improved the accuracy of the mean thickness estimates for synthetic cylinders and flat plates.

However, when you look at the image I attached to my comment above it seems like BinaryThinning is also kind of doing a crap job.. part of the joint space mask is totally missing in the local thickness map (top left corner) and that kind of issue must be what is throwing off the JSW stats.

I can try swapping it back to the old skeletonization algorithm and we can see if the JSW stats get better. It has its own problems when applied to perfectly flat or smooth shapes but maybe that's not a big deal with "real" data and maybe it will perform better in practice.

[njneeteson]

I switched back to skeletonize_3d and the results were even worse lol.. more of the mask was missing in the local thickness field than even before.

I found another library dedicated to skeletonization and it seems to be compatible with the other packages we are using. I'm going to try using it, I'll test a bit and probably push an update to my branch tomorrow.

[njneeteson]

OK... I have pushed an update to my branch. Now I don't use skeletonization at all. Instead of looping over voxels in the medial axis, I just loop over all voxels in the whole mask. The jited sphere-fitting function is fast so the computational optimization of using the medial axis isn't necessary and is introducing weird errors into the local thickness field.

Now that there is no skeletonization. I'm not really sure what could be wrong with the current local distance statistics function. Again maybe you will spot some errors and can let me know. But I've been looking at it on and off for a week and it seems to me like it is more correct than it was before.

Is it possible the actual problem is with the old method and the IPL method(s) have the same issue? Perhaps something to do with bad skeletonizations leading to inaccurate local thickness fields? When I look at the results now they seem pretty good:

For example for this JS mask, a JS.min of ~0.60 mm doesn't make sense to me unless the min_thickness parameter is set to 10*voxel_width. You can see on the left-hand side that the local thickness of those voxels goes all the way down to basically a single voxel "radius" which makes sense just from looking at the mask.

[mkuczyns]

I've taken look at the code and I can't see anything obviously wrong (at least right now, but I will take another look). From the image you've attached, wouldn't a JS.Min of ~0.6 mm make sense for the input mask? The minimum should be the smallest distance between bones (i.e., between the top and bottom edges of the mask) and not the smallest distance at the edges of the mask (if that makes sense). For this specific joint, the IPL derived JS.Min value is also ~0.60 mm.

I'll try downloading and running your branch to see how your updated DT method compares to known values from IPL. If the JS parameters are close to what I see from IPL, I'll work on merging your branch to main.

[njneeteson]

If you define JS.Min as the smallest distance between the two bones then probably 0.6 mm is reasonable. But when you apply Hildebrand's structure thickness algorithm to the JS mask it doesn't know that you only care about spheres that would fit in the mask and be touching both bones. It just fits spheres everywhere, including at those sharp corners, where the spheres get really small. So the minimum local thickness of the whole structure is not really the same as your definition of the minimum joint space width. If the values were matching up before I would guess that it was a coincidence due to overly sparse skeletonizations resulting in the local thickness field missing the corners of the JS mask where the local thickness gets small.

When you look at the report generated by IPL it seems like they are doing something to erode out the part of the JS mask that would stick out at the edges of the joint where the local thickness field would give you these "spurious" (for this application) smaller values. A potential fix would be to calculate the local thickness field of your whole JS mask then somehow do this erosion of the edge of the mask to get the smaller JS mask, then only take the local thickness values that are in this smaller region and then calculate mean, std, min, max, etc. on these values. Looking at the IPL report that seems like it has to be what is done in IPL otherwise those local thickness values in the contour plot for that mask don't make any sense to me.

I uploaded another notebook where I tried something like this but not quite the same:

JSW.mean : 1.558 mm
JSW.std  : 0.397 mm
JSW.min  : 0.738 mm
JSW.max  : 2.500 mm

With my attempt you get spurious larger values because big spheres fit into the mask outside of the joint. But you could follow a similar procedure with the JS mask and some kind of smaller eroded JS mask and probably get the same results as in IPL, maybe?

[njneeteson]

Tried my idea another way and got pretty much the same results:

JSW.mean : 1.558 mm
JSW.std  : 0.397 mm
JSW.min  : 0.738 mm
JSW.max  : 2.500 mm

So yeah, not really sure what is different between the IPL JSW algorithm and what we have in this repo that is giving different values. I definitely think that having smaller local thickness values in the corner of the JS mask makes sense and that something extra has to be done to get rid of them if they aren't present in the IPL results though.

[njneeteson]

Maybe we need to talk to Kathryn about how the local thickness function in IPL actually works, hah. It's possible they have some extra logic in there to handle voxels at the edges of a mask or in corners of a mask that is not in Hildebrand's paper and that I have not implemented in my version?

[mkuczyns]

I've also tried running your idea on a sample scan and the results are a bit closer to the IPL method than what we were previously running. One reason we might see a difference is that I am using the Python auto contour script but the IPL method's auto contour is not as "tight" to the bone. If I use one of the IPL generated masks, these are the results I get:

IPL JSW:

JS.W = 0.9698 mm JS.Min = 0.6070 mm JS.Max = 1.6389 mm

Your Updated JS Masking Method:

JS.W = 0.991 mm JS.Min = 0.543 mm JS.Max = 1.885 mm

Previous Python Method:

JS.W = 0.83576 mm JS.Min = 0.1214 mm JS.Max = 0.98625 mm

I'll have to dig a bit deeper to see why we see difference between your "new" method and the IPL method.

[njneeteson]

Nice well at least it is close. I'm going to do other stuff for now, let me know if after digging into it, there are any issues with my code that I need to fix. Once this is resolved and we're happy that the distance transform is working as expected, I will polish up the microarchitectural analysis code and then run it on a batch of images to see if that end is working properly.

[njneeteson]

I can also modify the current calc_stats or whatever function to optionally use this new idea if that would be useful for you. Would be pretty simple, just add an extra argument so you pass a mask to have the local thickness calculated on and then an optional second mask to use to mask over the local thickness field before calculating the stats. Then the workflow in a notebook or for a user is a lot cleaner.

[mkuczyns]

I've made some changes to the JS scripts in the main branch so I'll try to incorporate those with the new JS masking. That should make merging the branches a lot easier once we are happy with the DT.

I'll go through the IPL script for JS analysis today or tomorrow to make sure we're not missing something. I think we can setup a meeting with Kathryn in the next couple of weeks to show her what we have and if we're missing something in regards to the handling of edges in the JS. I'll send out an email to schedule something in the next few days.

[njneeteson]

OK sounds good. The changes needed were very minor so I modified the statistics calculating function to optionally take a sub-mask. Otherwise you would have to do a bunch of stuff in your JSW modules that would be kind of ugly (mask the local thickness yourself, set the minimum thickness yourself, etc.). And if anyone else uses this package maybe they'll find it useful for their own purposes.

I'm going to leave it alone for real now until it's time to merge or if fixes are needed.