Skip to content

Commit

Permalink
Change Evaluation Of Myelin Thickness And Area (#196)
Browse files Browse the repository at this point in the history 
* Change way to evaluate myelin thickness and area

Now, the thickness and area of the myelin takes is relative
to its corresponding axon.
Error handling for equivalent diameter and area values has
been added with the `_check_measures_are_relatively_valid`
method.

 Changes to be committed:
	modified:   AxonDeepSeg/morphometrics/compute_morphometrics.py

* Revert unit tests changes and update gratio formula

Gratio is axon_area / axonmyelin_area.

Another change is te remove raise ValueError when unexcepected measure are present.
Instead, a warning will be printed.

Changes to be committed:
	modified:   AxonDeepSeg/morphometrics/compute_morphometrics.py

* test_compute_mor[...]: Add UT For Bad Seg Case

* compute_mor[...]: Fix Warning Message For Measures

* test_compute_morp[...]: Change UT And Var Names

* [Bug] Fix definition of myelin thickness + test

* Fix test marker

* Fix Typo

* Clear Comments and Fix Warning Message
  • Loading branch information
@maf88 @mathieuboudreau
maf88 authored and mathieuboudreau committed Mar 30, 2019
1 parent 05cb434 commit 03d9492
Show file tree
Hide file tree
Showing 2 changed files with 154 additions and 53 deletions.
134 changes: 106 additions & 28 deletions AxonDeepSeg/morphometrics/compute_morphometrics.py
View file
@@ -1,6 +1,7 @@
# coding: utf-8

import os
from string import Template

# Scientific modules imports
import math
Expand Down Expand Up @@ -47,22 +48,24 @@ def get_axon_morphometrics(im_axon, path_folder, im_myelin=None):
:param im_myelin: Array: myelin binary mask, output of axondeepseg
:return: Array(dict): dictionaries containing morphometric results for each axon
"""
# TODO: externalize reading of pixel_size_in_micrometer.txt and input float
pixelsize = get_pixelsize(os.path.join(path_folder, 'pixel_size_in_micrometer.txt'))
stats_array = np.empty(0)
# Label each axon object
im_axon_label = measure.label(im_axon)
# Measure properties for each axon object
axon_objects = measure.regionprops(im_axon_label)

# Deal with myelin mask
if im_myelin is not None:
# sum axon and myelin masks

im_axonmyelin = im_axon + im_myelin
# Compute distance between each pixel and the background. Note: this distance is calculated from the im_axon,
# note from the im_axonmyelin image, because we know that each axon object is already isolated, therefore the
# distance metric will be more useful for the watershed algorithm below.

# Compute distance between each pixel and the background.
distance = ndi.distance_transform_edt(im_axon)
# local_maxi = feature.peak_local_max(distance, indices=False, footprint=np.ones((31, 31)), labels=axonmyelin)
# Note: this distance is calculated from the im_axon,
# note from the im_axonmyelin image, because we know that each axon
# object is already isolated, therefore the distance metric will be
# more useful for the watershed algorithm below.

# Get axon centroid as int (not float) to be used as index
ind_centroid = ([int(props.centroid[0]) for props in axon_objects],
Expand All @@ -74,20 +77,14 @@ def get_axon_morphometrics(im_axon, path_folder, im_myelin=None):
# Note: The value "i" corresponds to the label number of im_axon_label
im_centroid[ind_centroid[0][i], ind_centroid[1][i]] = i + 1

# markers = ndi.label(local_maxi)[0]
# Watershed segmentation of axonmyelin using distance map
im_axonmyelin_label = morphology.watershed(-distance, im_centroid, mask=im_axonmyelin)
# Measure properties of each axonmyelin object
axonmyelin_objects = measure.regionprops(im_axonmyelin_label)

# DEBUG
# from matplotlib import colors
# from matplotlib.pylab import *
# random_cmap = matplotlib.colors.ListedColormap(np.random.rand(256, 3))
# matshow(im_axon_label, fignum=1, cmap=random_cmap), show()
# import datetime
# savefig('fig_' + datetime.datetime.strftime(datetime.datetime.now(), '%Y%m%d%H%M%S%f') + '.png', format='png',
# transparent=True, dpi=100)
# Create list of the exiting labels
if im_myelin is not None:
axonmyelin_labels_list = [axm.label for axm in axonmyelin_objects]

# Loop across axon property and fill up dictionary with morphometrics of interest
for prop_axon in axon_objects:
Expand All @@ -111,30 +108,110 @@ def get_axon_morphometrics(im_axon, path_folder, im_myelin=None):
'solidity': solidity,
'eccentricity': eccentricity,
'orientation': orientation}

# Deal with myelin
if im_myelin is not None:
# Find label of axonmyelin corresponding to axon centroid
label_axonmyelin = im_axonmyelin_label[int(y0), int(x0)]
# TODO: use logger
# print(label_axonmyelin)
# print('x, y = {}, {}'.format(x0, y0))

if label_axonmyelin:
# Get corresponding index from axonmyelin_objects list
ind_axonmyelin = \
[axonmyelin_object.label for axonmyelin_object in axonmyelin_objects].index(label_axonmyelin)
myelin_diam = axonmyelin_objects[ind_axonmyelin].equivalent_diameter * pixelsize
myelin_area = axonmyelin_objects[ind_axonmyelin].area * (pixelsize ** 2)
stats['myelin_diam'] = myelin_diam
idx = axonmyelin_labels_list.index(label_axonmyelin)
prop_axonmyelin = axonmyelin_objects[idx]

_res1 = evaluate_myelin_thickness_in_px(prop_axon, prop_axonmyelin)
myelin_thickness = pixelsize * _res1

_res2 = evaluate_myelin_area_in_px(prop_axon, prop_axonmyelin)
myelin_area = (pixelsize ** 2) * _res2

axonmyelin_area = (pixelsize ** 2) * prop_axonmyelin.area

stats['myelin_thickness'] = myelin_thickness
stats['myelin_area'] = myelin_area
stats['gratio'] = np.sqrt(axon_area / myelin_area)
stats['axonmyelin_area'] = axonmyelin_area
stats['gratio'] = np.sqrt(axon_area / axonmyelin_area)
else:
# TODO: use logger
print('WARNING: Myelin object not found for axon centroid [{},{}]'.format(y0, x0))
print(
"WARNING: Myelin object not found for axon" +
"centroid [y:{0}, x:{1}]".format(y0, x0)
)

stats_array = np.append(stats_array, [stats], axis=0)

return stats_array

def evaluate_myelin_thickness_in_px(axon_object, axonmyelin_object):
"""
Returns the equivalent thickness of a myelin ring around an axon of a
given equivalent diameter (see note [1] below). The result is in pixels.
:param axon_object (skimage.measure._regionprops): object returned after
measuring a axon labeled region
:param axonmyelin_object (skimage.measure._regionprops): object returned after
measuring a axon with myelin labeled region
[1] According to https://scikit-image.org/docs/dev/api/skimage.measure.html?highlight=region%20properties#regionprops,
the equivalent diameter is the diameter of a circle with the same area as
the region.
"""
warn_if_measures_are_unexpected(
axon_object,
axonmyelin_object,
"equivalent_diameter"
)

axon_diam = axon_object.equivalent_diameter
axonmyelin_diam = axonmyelin_object.equivalent_diameter
return (axonmyelin_diam - axon_diam)/2

def evaluate_myelin_area_in_px(axon_object, axonmyelin_object):
"""
Returns the myenlinated axon area minus the axon area.
:param axon_object (skimage.measure._regionprops): object returned after
measuring an axon labeled region
:param axonmyelin_object (skimage.measure._regionprops): object returned after
measuring a axon with myelin labeled region
"""
warn_if_measures_are_unexpected(
axon_object,
axonmyelin_object,
"area"
)
return axonmyelin_object.area - axon_object.area

def warn_if_measures_are_unexpected(axon_object, axonmyelin_object, attribute):
"""
Calls the `_check_measures_are_relatively_valid` function and if return
value is False, print a warning.
"""
checked = _check_measures_are_relatively_valid(axon_object, axonmyelin_object, attribute)
if checked is False:
x_a, y_a = axon_object.centroid
data = {
"attribute": attribute,
"axon_label": axon_object.label,
"x_ax": x_a,
"y_ax": y_a,
"axonmyelin_label": axonmyelin_object.label,
}

warning_msg = Template(
"Warning, axon #$axon_label at [y:$y_ax, x:$x_ax] and " +
"corresponding myelinated axon #$axonmyelin_label " +
"have unexpected measure values for $attribute attributest."
)
print(warning_msg.safe_substitute(data))

def _check_measures_are_relatively_valid(axon_object, axonmyelin_object, attribute):
"""
Checks if the attribute is positive and if the myelinated axon has a greater value
"""
val_axon = getattr(axon_object, attribute)
val_axonmyelin = getattr(axonmyelin_object, attribute)
if val_axon > 0 and val_axonmyelin > 0 and val_axonmyelin > val_axon:
return True
else:
return False

def save_axon_morphometrics(path_folder, stats_array):
"""
Expand Down Expand Up @@ -169,7 +246,7 @@ def load_axon_morphometrics(path_folder):
def draw_axon_diameter(img, path_prediction, pred_axon, pred_myelin):
"""
:param img: sample grayscale image (png)
:param path_prediction: full path to the segmented file (*_seg-axonmyelin.png)
:param path_prediction: full path to the segmented file (*_seg-axonmyelin.png)
from axondeepseg segmentation output
:param pred_axon: axon mask from axondeepseg segmentation output
:param pred_myelin: myelin mask from axondeepseg segmentation output
Expand Down Expand Up @@ -273,3 +350,4 @@ def write_aggregate_morphometrics(path_folder, aggregate_metrics):
"directory \"{1}\".\n".format('aggregate_morphometrics.txt', path_folder)))

raise

73 changes: 48 additions & 25 deletions test/morphometrics/test_compute_morphometrics.py
View file
Expand Up @@ -105,7 +105,7 @@ def test_get_axon_morphometrics_with_myelin_mask(self):
flatten=True)

stats_array = get_axon_morphometrics(pred_axon, path_folder, im_myelin=pred_myelin)

assert stats_array[1]['gratio'] == pytest.approx(0.74, rel=0.01)

@pytest.mark.unit
Expand All @@ -115,30 +115,32 @@ def test_get_axon_morphometrics_with_myelin_mask_simulated_axons(self):
'__test_files__',
'__test_simulated_axons__',
'SimulatedAxons.png')

gratio_sim = [
0.9,
0.8,
0.7,
0.6,
0.5,
0.4,
0.3,
0.2,
0.1
]

axon_diam_sim = [
100,
90,
80,
70,
60,
46,
36,
24,
12
]

gratio_sim = np.array([
0.9,
0.8,
0.7,
0.6,
0.5,
0.4,
0.3,
0.2,
0.1
])

axon_diam_sim = np.array([
100,
90,
80,
70,
60,
46,
36,
24,
12
])

myelin_thickness_sim = (axon_diam_sim/2)*(1/gratio_sim-1)

# Read paths and compute axon/myelin masks
pred = scipy_imread(path_pred)
Expand All @@ -152,7 +154,28 @@ def test_get_axon_morphometrics_with_myelin_mask_simulated_axons(self):
for ii in range(0,9):
assert stats_array[ii]['gratio'] == pytest.approx(gratio_sim[ii], rel=0.1)
assert stats_array[ii]['axon_diam'] == pytest.approx(axon_diam_sim[ii], rel=0.1)
assert stats_array[ii]['myelin_thickness'] == pytest.approx(myelin_thickness_sim[ii], rel=0.1)

@pytest.mark.unit
def test_get_axon_morphometrics_with_unexpected_myelin_mask_simulated_axons(self):
path_pred = os.path.join(
self.testPath,
'__test_files__',
'__test_simulated_axons__',
'SimulatedAxons.png')

# Read paths and compute axon/myelin masks
pred = scipy_imread(path_pred, flatten=True)
pred_axon = pred > 200
unexpected_pred_myelin = np.zeros(pred.shape)
path_folder, file_name = os.path.split(path_pred)

# Compute axon morphometrics
stats_array = get_axon_morphometrics(pred_axon,path_folder,im_myelin=unexpected_pred_myelin)
for axon_prop in stats_array:
assert axon_prop['myelin_thickness'] == pytest.approx(0.0, rel=0.01)
assert axon_prop['myelin_area'] == pytest.approx(0.0, rel=0.01)
assert axon_prop['gratio'] == pytest.approx(1.0, rel=0.01)

# --------------save and load _axon_morphometrics tests-------------- #
@pytest.mark.unit
Expand Down

0 comments on commit 03d9492

Please sign in to comment.