Example of watershed in slam

# Authors:
# Lucile Hashimoto lucile-hashimoto
# Guillaume Auzias <guillaume.auzias@univ-amu.fr>

# License: MIT
# sphinx_gallery_thumbnail_number = 2

importation of slam modules

import os
import numpy as np
import slam.texture as stex
import slam.io as sio
import slam.watershed as swat

loading an examplar mesh and corresponding texture

path_to_mesh = "../examples/data/example_mesh.gii"
path_to_mask = None
path_to_output = None
mesh = sio.load_mesh(path_to_mesh)

set the threshold values for the watershed

thresh_dist_perc, thresh_ridge, thresh_area_perc = 10, 0.0001, 3

compute curvature, dpf and voronoi

mean_curvature, dpf_star, voronoi = swat.compute_mesh_features(mesh)

        Computing the curvature

Calculating vertex normals .... Please wait
Finished calculating vertex normals
Calculating curvature tensors ... Please wait
Finished Calculating curvature tensors
Calculating Principal Components ... Please wait
Finished Calculating principal components

        Computing the DPF

  Computing Laplacian
    Computing mesh weights of type conformal
    -edge length threshold needed for  0  values =  0.0  %
    -number of Nan in weights:  0  =  0.0  %
    -number of Negative values in weights:  936  =  6.706792777300086  %
    -nb Nan in Laplacian :  0
    -nb Inf in Laplacian :  0

        Computing Voronoi's vertex

    -percent polygon with obtuse angle  29.4067067927773

define the exclusion mask (cingular pole)

if path_to_mask is not None:
    mask = sio.load_texture(path_to_mask).darray[0]
else:
    mask = None

extract sulcal pits and associated basins

basins, ridges, adjacency = swat.watershed(
    mesh, voronoi, dpf_star, thresh_dist_perc, thresh_ridge, thresh_area_perc, mask)

Computing watershed by flooding...
Thresholds provided:
-Distance between 2 pits: 10 % <=> 3.3720850389113624 mm
-Ridge height: 0.0001
-Basin area: 3 % <=> 193.10651237927843 mm2

print basic info

print("nb basins:", len(basins))
print("nb ridges:", len(ridges))

nb basins: 15
nb ridges: 39

print statistics regarding the area of basins

bas_area = swat.get_basins_attribute(basins, attribute='basin_area')
print('average basins area:', np.mean(bas_area))
print('range of basins area:', min(bas_area), '-', max(bas_area))

average basins area: 429.12558306506315
range of basins area: 194.51144922010275 - 1145.6073016360745

print statistics regarding the height of ridges

rid_height = swat.get_ridges_attribute(ridges, attribute='ridge_depth_diff_min')
print('average ridges height:', np.mean(rid_height))
print('range of ridges height:', min(rid_height), '-', max(rid_height))

average ridges height: 0.0030232953255292143
range of ridges height: 0.0 - 0.012923161294396496

generate the textures from watershed outputs

atex_labels, atex_pits, atex_ridges = swat.get_textures_from_dict(
    mesh, basins, ridges)

print the labels present in the texture of basin labels

print(np.unique(atex_labels))

[ 0  1  3  4  6  8 10 12 13 15 16 18 20 21 26]

generate the texture of the boundaries between basins from watershed outputs

atex_boundaries = swat.get_texture_boundaries_from_dict(mesh, ridges)

if path_to_output is not None:
    # texture of curvature
    curv_tex = stex.TextureND(darray=mean_curvature)
    sio.write_texture(curv_tex, os.path.join(path_to_output, "mean_curvature.gii"))
    # texture of depth
    dpf_tex = stex.TextureND(darray=dpf_star)
    sio.write_texture(dpf_tex, os.path.join(path_to_output, "dpf_star.gii"))
    # texture of voronoi
    voronoi_tex = stex.TextureND(darray=voronoi)
    sio.write_texture(voronoi_tex, os.path.join(path_to_output, "voronoi.gii"))
    # texture of labels
    tex_labels = stex.TextureND(darray=atex_labels)
    sio.write_texture(tex_labels, os.path.join(path_to_output, "basins.gii"))
    # texture of pits
    tex_pits = stex.TextureND(darray=atex_pits)
    sio.write_texture(tex_pits, os.path.join(path_to_output, "pits.gii"))
    # texture of ridges
    tex_ridges = stex.TextureND(darray=atex_ridges)
    sio.write_texture(tex_ridges, os.path.join(path_to_output, "ridges.gii"))
    # texture of ridges vertices
    texture_boundaries = stex.TextureND(darray=atex_boundaries)
    sio.write_texture(tex_labels, os.path.join(path_to_output, "bondaries.gii"))

VISUALIZATION USING plotly

import slam.plot as splt
# set the parameter to shift the pits and ridges outwards for better visualization
out_shift=0.2
# create a texture combining the DPF* with basins boundaries
tex_plot = dpf_star.copy()
min_dpf = min(dpf_star)
tex_plot[atex_boundaries==1] = min_dpf

display_settings = {}
display_settings['colorbar_label'] = 'DPF*'
mesh_data = {}
mesh_data['vertices'] = mesh.vertices
mesh_data['faces'] = mesh.faces
intensity_data = {}
intensity_data['values'] = tex_plot
intensity_data["mode"] = "vertex"
fig = splt.plot_mesh(
    mesh_data=mesh_data,
    intensity_data=intensity_data,
    display_settings=display_settings)
# add the ridges to the plot
ridges_coords = splt.coords_outward_shift(atex_ridges==1, mesh, out_shift=0.5)
trace_ridges = splt.plot_points(
    ridges_coords,
    marker={"size": 6, "color": "white", "line":{"color":"black", "width":2}},
)
fig.add_trace(trace_ridges)
# add the pits to the plot
pits_coords = splt.coords_outward_shift(atex_pits==1, mesh, out_shift=1)
trace_pits = splt.plot_points(
    pits_coords,
    marker={"size": 6, "color": "red", "line":{"color":"black", "width":2}}
)
fig.add_trace(trace_pits)
fig.show()
fig

# create a texture showing the basins label and boundaries
tex_plot = atex_labels.copy()
tex_plot[atex_boundaries==1] = -1
print(np.unique(tex_plot))

display_settings = {}
display_settings['colorbar_label'] = 'Basins label'
display_settings["colorscale"]= "BrBg"
intensity_data = {}
intensity_data['values'] = tex_plot
intensity_data["mode"] = "vertex"
fig2 = splt.plot_mesh(
    mesh_data=mesh_data,
    intensity_data=intensity_data,
    display_settings=display_settings)
# add the ridges to the plot
fig2.add_trace(trace_ridges)
# add the pits to the plot
fig2.add_trace(trace_pits)
fig2.show()
fig2

[-1  0  1  3  4  6  8 10 12 13 15 16 18 20 21 26]

extract corresponding sulcal graph

import slam.sulcal_graph as ssg
g = ssg.get_sulcal_graph(adjacency, basins, ridges)
g = ssg.add_coords_attribute(g, mesh,
                                    attribute_vert_index='pit_index',
                                    new_attribute_key='3dcoords')

mesh_data["opacity"]=0.5
fig3 = splt.plot_mesh(
    mesh_data=mesh_data,
    intensity_data=intensity_data,
    display_settings=display_settings)
fig3 = splt.plot_graph_on_mesh(g, mesh, vertex_index_attribute='pit_index', out_shift=1, show_edge_ridge=True, fig=fig3)
# add the ridges to the plot
fig3.add_trace(trace_ridges)
fig3.show()
fig3

{'pit_index': 530, 'basin_vertices': [530, 535, 2172, 529, 524, 2131, 2173, 531, 947, 964, 941, 911, 2199, 2174, 559, 912, 2132, 2187, 2198, 910, 909, 2158, 546, 2091, 561, 513, 2074, 2200, 871, 948, 537, 927, 860, 2211, 2133, 2188, 894, 861, 895, 862, 2216, 863, 875, 846, 966, 2204, 2212, 2104, 496, 873, 514, 913, 2038, 560, 2218, 847, 2134, 819, 876, 965, 2135, 2159, 2213, 516, 967, 949, 462, 874, 503, 548, 794, 2161, 2138, 2215, 2136, 2205, 463, 497, 968, 806, 782, 929, 884, 2163, 2219, 356, 2217, 831, 486, 464, 2220, 2227, 450, 155, 553, 2214, 2230, 848, 440, 952, 795, 456, 2177, 2190, 878, 741, 2232, 2231, 781, 1925, 774, 547, 896, 2242, 747, 2191, 808, 2221, 914, 942, 566, 885, 2250, 536, 2264, 2235, 769, 797, 850, 915, 2251, 552, 578, 515, 2265, 743, 879, 851, 2233, 886, 2266, 621, 770, 596, 2234, 587, 880, 488, 605, 562, 615, 2222, 579, 809, 852, 629, 853, 604, 2244, 640, 489, 748, 2252, 2236, 641, 458, 2223, 2243, 639, 608, 2267, 2268, 427, 2280, 403, 411, 393, 382, 370, 1529, 2092, 518, 2105, 517, 2139, 504, 2137, 2149, 877, 2160, 2162, 2176, 2178, 2179, 2193], 'pit_depth': np.float64(-0.013045193486841064), 'basin_area': np.float64(483.1587592928476), 'basin_label': 0, '3dcoords': array([ 6.71122503, -2.44120479,  0.5532046 ])}