test-2-and-3-D.py

# Copyright CNRS/Inria/UNS
# Contributor(s): Eric Debreuve (since 2018)
#
# eric.debreuve@cnrs.fr
#
# This software is governed by the CeCILL  license under French law and
# abiding by the rules of distribution of free software.  You can  use,
# modify and/ or redistribute the software under the terms of the CeCILL
# license as circulated by CEA, CNRS and INRIA at the following URL
# "http://www.cecill.info".
#
# As a counterpart to the access to the source code and  rights to copy,
# modify and redistribute granted by the license, users are provided only
# with a limited warranty  and the software's author,  the holder of the
# economic rights,  and the successive licensors  have only  limited
# liability.
#
# In this respect, the user's attention is drawn to the risks associated
# with loading,  using,  modifying and/or developing or reproducing the
# software by the user in light of its specific status of free software,
# that may mean  that it is complicated to manipulate,  and  that  also
# therefore means  that it is reserved for developers  and  experienced
# professionals having in-depth computer knowledge. Users are therefore
# encouraged to load and test the software's suitability as regards their
# requirements in conditions enabling the security of their systems and/or
# data to be ensured and,  more generally, to use and operate it in the
# same conditions as regards security.
#
# The fact that you are presently reading this means that you have had
# knowledge of the CeCILL license and that you accept its terms.

import matplotlib as mpl_

mpl_.use('TkAgg')

from sklgraph.skl_fgraph import skl_graph_t
from sklgraph.skl_graph import plot_mode_e
from sklgraph.skl_map import skl_map_t

import pathlib as ph_
import sys as sy_
import tempfile as tf_

import matplotlib.pyplot as pl_
import numpy as np_
import scipy.ndimage as im_
import skimage.measure as ms_

# import sknw as sk_


# --- Set main parameters
if sy_.argv.__len__() == 1:
    print("Call with arguments: horse|bug|2d|3d [noplot]")
    sy_.exit(0)

data = sy_.argv[1]
print(data)
test_mode = (sy_.argv.__len__() > 2) and (sy_.argv[2] == "noplot")
space_dim = None

# --- Generation of a (random) object map
if data == "horse":
    import skimage.data as sd_
    import skimage.util as su_

    obj_map = su_.invert(sd_.horse())
    space_dim = 2
#
elif data == "bug":
    bug_path = ph_.Path(__file__).parent
    bug_path = bug_path / ".." / "dev" / "__usecases__" / "obj_map.npz"
    if bug_path.is_file():
        obj_map = np_.load(bug_path)["obj_map"]
        space_dim = obj_map.ndim
    else:
        raise FileNotFoundError(bug_path)
#
else:
    if data == "2d":
        space_dim = 2
        array_shape = (100, 100)
        threshold_fct = lambda img: img.median()
        struct_elm_shape = (3, 3)
    elif data == "3d":
        space_dim = 3
        array_shape = (80, 80, 50)
        threshold_fct = lambda img: 0.75 * image.max()
        struct_elm_shape = (3, 3, 3)
    else:
        raise ValueError(f"{data}: Unknown data specification")

    tmp_folder = None
    with tf_.NamedTemporaryFile() as doc_accessor:
        doc_name = ph_.Path(doc_accessor.name)
        tmp_folder = doc_name.parents[0]
    if tmp_folder is not None:
        seed_lock = tmp_folder.joinpath("skl_graph_demo_seed.lock")
        if not seed_lock.exists():
            np_.random.seed(0)
            seed_lock.touch()

    image = np_.random.random(array_shape)
    image = im_.gaussian_filter(image, 1)

    multi_map = image > threshold_fct(image)
    lbl_multi_map, _ = im_.label(
        multi_map,
        structure=np_.ones(struct_elm_shape, dtype=np_.uint8),
        output=np_.uint,
    )

    cc_props = ms_.regionprops(lbl_multi_map)
    largest_cc_idx = np_.argmax([cc_prop.area for cc_prop in cc_props]).item()
    lbl_multi_map[lbl_multi_map != cc_props[largest_cc_idx].label] = 0
    obj_map = lbl_multi_map > 0

# --- From object map to skeleton graph
np_.savez_compressed(".\\obj_map", obj_map=obj_map)
skeleton = skl_map_t.FromShapeMap(obj_map, store_widths=True)
skel_graph = skl_graph_t.FromSkeleton(skeleton)

# sknw_graph = sk_.build_sknw(skeleton.map)
# print(dir(sknw_graph))

# --- Some info about the skeleton graph
print(
    f"Obj map area={np_.count_nonzero(obj_map)}\n\n"
    f"Validity={skel_graph.is_valid}\n\n"
    f"N nodes={skel_graph.n_nodes}\n"
    f"N edges={skel_graph.n_edges}\n"
    f"Highest degree={skel_graph.highest_degree}/{skel_graph.highest_degree_w_nodes}\n\n"
    f"Length={skel_graph.length}<-{skel_graph.edge_lengths}\n"
    f"Width=Hom.{skel_graph.reduced_width()}/Het.{skel_graph.heterogeneous_reduced_width()}<-{skel_graph.edge_reduced_widths()}\n"
    f"Area as WxL={skel_graph.reduced_width() * skel_graph.length}\n"
    f"Area as WW Length={skel_graph.ww_length}<-{skel_graph.edge_ww_lengths}\n"
)

# --- No display is testing
if test_mode:
    sy_.exit(0)

# --- Checking correctness of skeleton graph via various rebuilt maps
rebuilt_skel = skel_graph.RebuiltSkeletonMap()
part_map_for_check = skeleton.PartMap()
# Keep next line before its next one
part_map_for_check[part_map_for_check == skeleton.invalid_n_neighbors] = 0
part_map_for_check[part_map_for_check > 3] = 3

print(
    f"Rebuilt skeleton = Original one? Part map? "
    f"{np_.array_equal(rebuilt_skel > 0, skeleton.map)} "
    f"{np_.array_equal(rebuilt_skel, part_map_for_check)}"
)

if space_dim == 3:
    rebuilt_skel = rebuilt_skel.max(axis=2)

# --- Checking correctness of skeleton graph via rebuilt object map
obj_map = obj_map.astype(np_.uint8)
rebuilt_map = skel_graph.RebuiltObjectMap()
if space_dim == 3:
    obj_map = obj_map.max(axis=2)
    rebuilt_map = rebuilt_map.max(axis=2)

# --- Show various images and plot the skeleton graph
___, axes = pl_.subplots(1, 3)
axes[0].matshow(obj_map)
axes[1].matshow(rebuilt_map)
axes[2].matshow(obj_map - rebuilt_map)
axes[0].set_title("Object")
axes[1].set_title("Rebuilt")
axes[2].set_title("Object - Rebuilt")
for axis in axes:
    axis.set_axis_off()
pl_.tight_layout(pad=0.3)

___, axes = pl_.subplots(1, 3)
axes[0].matshow(rebuilt_skel + obj_map)
axes[1].matshow(rebuilt_skel + rebuilt_map)
axes[2].matshow(rebuilt_skel)
axes[0].set_title("Rebuilt Skel Over Obj")
axes[1].set_title("Rebuilt Skel Over Rebuilt Obj")
axes[2].set_title("Rebuilt Skel")
for axis in axes:
    axis.set_axis_off()
pl_.tight_layout(pad=0.3)

if space_dim == 2:
    skel_graph.Plot(mode=plot_mode_e.Networkx, should_block=False)
    skel_graph.Plot(mode=plot_mode_e.Graphviz, should_block=False)
skel_graph.Plot(mode=plot_mode_e.SKL_Curve, w_directions=True, should_block=False)
skel_graph.Plot(should_block=False)

# From: https://github.com/Image-Py/sknw (previously: https://github.com/yxdragon/sknw)
# sk_.draw_graph(obj_map, sknw_graph, cn=255, ce=128)
# pl_.matshow(obj_map, cmap="gray")
# pl_.title("SKNW")
# pl_.gca().set_axis_off()
# pl_.tight_layout(pad=0.3)

pl_.show()