Source code for soma.aimsalgo.meshwatershedtools

#! /usr/bin/env python
############################################################################
# This software and supporting documentation are distributed by
# CEA/NeuroSpin, Batiment 145, 91191 Gif-sur-Yvette cedex, France.
# This software is governed by the CeCILL license version 2 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 version 2 as circulated by CEA, CNRS
# and INRIA at the following URL "http://www.cecill.info".
############################################################################

"""
This modles features the following:

* generate a reduced profile by basins using the watershed algorithm
* merge the items if necessary (measured criteria validating the basins)

Main dependencies: PyAims library
"""


#----------------------------Imports-------------------------------------------


# import system
from __future__ import absolute_import
from __future__ import print_function
import numpy
import sys

# soma import
from soma import aims, aimsalgo
from six.moves import range

#----------------------------Functions-----------------------------------------


def _clean_after_merge(labels, parents, valid):
    """Remove the non-valid components and reorder the ROI list.append

    Parameters
    ----------
    labels:
        labels of basins of watershed (after merge)
    parents:
        parent local according to local maximum (after merge)
    valid:
        measured criteria validating the basins

    Returns
    -------
    labels:
        clean labels of basins of watershed
    parents:
        clean parent local according to local maximum
    """
    k = numpy.size(valid)
    if k != 0:
        ## remove invalid null components
        for j in range(k):
            if not valid[parents[j]]:
                parents[j] = j

        iconvert = numpy.nonzero(valid)
        iconvert = numpy.reshape(iconvert, numpy.size(iconvert))
        convert = -numpy.ones(k).astype('i')
        aux = (numpy.cumsum(valid.astype('i')) - 1).astype('i')
        convert[valid] = aux[valid]
        k = numpy.size(iconvert)

        q = numpy.nonzero(labels > -1)
        q = numpy.reshape(q, numpy.size(q))
        labels[q] = convert[labels[q]]
        parents = convert[parents[iconvert]]
        parents = parents[:k]

        return labels, parents


def _merge_ascending(labels, parents, valid):
    """Remove the non-valid items by including them in their parents when
    it exists.

    Parameters
    ----------
    labels:
        labels of basins of watershed
    parents:
        parent local according to local maximum
    valid:
        measured criteria validating the basins

    Returns
    -------
    labels:
        clean labels of basins of watershed
    parents:
        clean parents according to local maximum
    """
    labels = numpy.copy(labels)
    parents = numpy.copy(parents)
    for j in range(len(valid)):
        if valid[j] == 0:
            fj = parents[j]
            if fj != j:
                parents[parents == j] = fj
                labels[labels == j] = fj
            else:
                valid[j] = 1

    labels, parents = _clean_after_merge(labels, parents, valid)

    return labels, parents


[docs]def watershedWithBasinsMerging( tex, mesh, size_min=0, depth_min=0, tex_threshold=0.01, mode="and"): """Generate a texture of merged basins: watershed texture. The basins are merged according to two criteria: (1) the size of basins (2) the depth of basins Parameters ---------- tex: (TimeTexture_S16) texture of boundaries between regions mesh: (aimsTimeSurface_3) associated mesh size_min: (int) number of basins > size_min depth_min: (float) number of basins > depth_min tex_threshold: (float) threshold on the input intensity texture mode: (str) two cases: (1) and --> merge basins with its parent if size < k_size and depth < k_depth (2) or --> merge basins with its parent if size < k_size or depth < k_depth Returns ------- output_tex: (TimeTexture_S16) watershed results according to the thresholds indicated by k_size and k_depth """ ## Watershed: idxWt, depthWt, majorWt, labelWt = aimsalgo.meshWatershed( mesh, tex, tex_threshold) # idxWt : array of shape (nbasins) indices of the vertices that are # local maxima # depthWt: array of shape (nbasins) topological the depth of the basins # depth[idx[i]] = q means that idx[i] is a q-order maximum # This is also the diameter of the basins associated # with local maxima # majorWt: array of shape (nbasins) # label of the maximum which dominates each local maximum # i.e. it describes the hierarchy of the local maxima # labelWt: array of shape (number of vertices) # labelling of the vertices according to their basin data = tex[0].arraydata() idxW = idxWt[0].data().arraydata() majorW = majorWt[0].data().arraydata() labelW = labelWt[0].data().arraydata() if size_min != 0 or depth_min != 0: size = numpy.asarray([len(numpy.where(labelW == x)[0]) for x in range(idxW.size)]) print("criteria size: ", size) # Blobs blobindices = aimsalgo.blobsHeights( mesh, tex[0].data(), labelWt[0].data()) blobheights = data[blobindices.arraydata()] blobdepths = data[idxW] - blobheights print("criteria depth: ", blobdepths) # criteria to regroup basins valid_size = size > size_min valid_depth = blobdepths > depth_min if mode == "or": valid = valid_size * valid_depth elif mode == "and": valid = valid_size + valid_depth else: raise ValueError('invalid mode') # Remove the non-valid items by including them in their parents when # it exists. labelW, majorW = _merge_ascending(labelW, majorW, valid) # create the final texture with the items valid output_tex = aims.TimeTexture_S16() output_text = output_tex[0] output_text.assign(labelW + 1) return output_tex