Source code for soma.aims.colormaphints

# -*- coding: utf-8 -*-
#  This software and supporting documentation are distributed by
#      Institut Federatif de Recherche 49
#      CEA/NeuroSpin, Batiment 145,
#      91191 Gif-sur-Yvette cedex
#      France
#
# This software is governed by the CeCILL-B 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-B 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-B license and that you accept its terms.

from __future__ import absolute_import
from __future__ import print_function
from soma import aims
import numpy
import types
import sys
from six.moves import range
try:
    from soma import aimsalgo
except:
    aimsalgo = None


[docs]def checkVolume(vol): '''Checks colormap-related clues in a volume, and tries to determine whether it is an anatomical volume, a diffusion volume, a functional volume, or a labels volume. This is determined as "likelihoods" for each class (based on a pure empirical heurisrtic), based on, mainly, the histogram, voxel type, and voxel sizes. ''' def _histogram(vol, bins): histo = getattr(aims, 'RegularBinnedHistogram_' + dtype) if type(bins) is not int: m, M = bins[0], bins[-1] if dtype in ('U8', 'S8'): m = chr(m) M = chr(M) bins = len(bins) - 1 hg = histo(bins) hg.doit(vol, m, M) else: hg = histo(bins) hg.doit(vol) d = hg.data() h = numpy.array(d.volume(), copy=False).reshape(d.dimX()) if dtype in ('U8', 'S8'): step = float(ord(hg.maxDataValue()) - ord(hg.minDataValue())) / hg.bins() his = (h, numpy.arange(ord(hg.minDataValue()), ord(hg.maxDataValue()) + step, step)) else: step = float(hg.maxDataValue() - hg.minDataValue()) / hg.bins() his = (h, numpy.arange(hg.minDataValue(), hg.maxDataValue() + step, step)) return his def _unique(vol): hg = getattr(aims, 'RegularBinnedHistogram_' + dtype)() # abort if more than 100000 values (this will not be labels) return hg.unique(vol, 100000).arraydata() if aimsalgo is not None: # aims histogram is 15 times faster than numpy... histogram = _histogram # and unique is only 5 times faster than numpy unique = _unique else: histogram = numpy.histogram unique = lambda vol: numpy.unique(numpy.array(vol, copy=False)) hints = {} # check data type tname = vol.__class__.__name__ if tname.startswith('Volume_'): dtype = tname[7:] elif tname.startswith('AimsData_'): dtype = tname[9:] vol = vol.volume() else: raise ValueError('input parameter should be a volume') hints['data_type'] = dtype if dtype not in ('U8', 'S8', 'U16', 'S16', 'U32', 'S32', 'U64', 'S64', 'FLOAT', 'DOUBLE'): # not a scalar, we cannot perform scalar things, and the volume will have # no colormap, so we are done return hints try: u = unique(vol) maxv = u[-1] minv = u[0] except: maxv = vol.max() minv = vol.min() u = None hints['binary_volume'] = False hints['int_type'] = False hints['labels_likelihood'] = 0. if u is not None: if len(u) == 2: hints['binary_volume'] = True else: hints['binary_volume'] = False if dtype[0] == 'S' or dtype[0] == 'U' and dtype[1:] in ('8', '16', '32', '64'): hints['int_type'] = True else: if len(numpy.where(u.astype('int32') == u)) != len(u): hints['int_type'] = False hints['labels_likelihood'] = 0. else: hints['int_type'] = True if hints['int_type']: labell = maxv / 250. - 1. labell = numpy.exp(-labell) / (numpy.exp(labell) + numpy.exp(-labell)) hints['labels_likelihood'] = labell if minv < 0: hints['negative_values'] = True else: hints['negative_values'] = False # determine background color (most frequent value) if u is not None: nbins = min(len(u), 200) else: nbins = 200 his = histogram(vol, nbins) maxhis = numpy.argmax(his[0]) if u is not None: subval = numpy.where((u < his[1][maxhis + 1]).__and__( u >= his[1][maxhis]))[0] if len(subval) == 0: subval = [0] # FIXME: take the next lower used bin if len(subval) == 1: hints['background'] = u[subval[0]] else: # more precise histogram bins = numpy.hstack((u[subval], u[subval[-1] + 1])) h = histogram(vol, bins) hints['background'] = h[1][numpy.argmax(h[0])] maxhis = numpy.where(his[1] <= hints['background'])[0][-1] else: # float values hints['background'] = minv # determine cmap bounds # "positive" values chis = numpy.cumsum(his[0][maxhis:]) chis = chis / float(chis[-1]) # proportional cummulated histo maxcmap = numpy.where(chis < 0.99)[0] if len(maxcmap) >= 1: maxcmap = maxcmap[-1] if chis[maxcmap] < 0.98: maxcmap += 1 # add the next bin if it is too big else: maxcmap = len(chis) if maxcmap < len(chis): hints['histo_99percent_positive'] = his[1][maxhis + maxcmap] else: hints['histo_99percent_positive'] = maxv # correct the labels likelihood hints['labels_likelihood'] *= float(maxcmap) / len(chis) if hints['binary_volume']: hints['labels_likelihood'] = 1. # "negative" values chis = numpy.cumsum(his[0][:maxhis + 2][-1:0:-1]) chis = chis / float(chis[-1]) # proportional cummulated histo maxcmap = numpy.where(chis < 0.99)[0] if len(maxcmap) >= 1: maxcmap = maxcmap[-1] if chis[maxcmap] < 0.98: maxcmap += 1 # add the next bin if it is too big else: maxcmap = len(chis) if maxcmap < len(chis): hints['histo_99percent_negative'] = his[1][maxhis - maxcmap] else: hints['histo_99percent_negative'] = minv hdr = vol.header() dims = [vol.getSizeX(), vol.getSizeY(), vol.getSizeZ(), vol.getSizeT()] vs = hdr['voxel_size'] # determine if it could be an anatomical / functional volume vsmin = numpy.min(vs[:3]) vsmax = numpy.max(vs[:3]) def anatlikelihoodcurve(vsize): if vsize < 1.: return 1. elif vsize > 4.: return 0.2 else: return max(1. - (vsize - 1.) / 3., 0.2) def difflikelihoodcurve(vsize): if vsize < 1.: return 0.2 elif vsize > 6.: return 0.1 elif vsize < 2.: return max(vsize - 1., 0.2) else: return max(1. - (vsize - 2.) / 4., 0.1) def funclikelihoodcurve(vsize): if vsize < 1.: return 0.1 elif vsize > 3.: return 1. else: return max((vsize - 1) / 2., 0.1) def normalize(vec): return vec / numpy.sqrt((vec ** 2).sum()) likes = normalize(numpy.array([anatlikelihoodcurve(vsmin), difflikelihoodcurve(vsmin), funclikelihoodcurve(vsmin), 0.])) likes *= [anatlikelihoodcurve(vsmax), difflikelihoodcurve(vsmax), funclikelihoodcurve(vsmax), 0.] likes = normalize(likes) likes[3] = hints['labels_likelihood'] if not hints['int_type']: if hints['negative_values']: likes *= [0.2, 0.3, 1., 0.] else: if dims[3] > 1: likes *= [0.1, 0.5, 2., 0.] else: likes *= [0.3, 0.3, 2., 0.] else: ll = 1. - hints['labels_likelihood'] if hints['negative_values']: if dims[3] > 1: likes *= [0.2 * ll, 0.5 * ll, ll, 1.] else: likes *= [0.5 * ll, 0.2 * ll, ll, 1.] else: likes *= [ll, ll, ll, 1.] likes = normalize(likes) hints['volume_contents_likelihoods'] = list(likes) hints['anatomical_likelihood'] = likes[0] hints['diffusion_likelihood'] = likes[1] hints['functional_likelihood'] = likes[2] hints['labels_likelihood'] = likes[3] # get specified palette information if any try: print('hdr:', type(hdr)) print(hdr) pal = hdr['palette'] p = pal['palette'] hints['palette'] = p except: # KeyError: # FIXME: with python3 we can get an AttribteError pass return hints
anatomicalColormaps = [('B-W LINEAR', (1., 1., 1.)), ('Blue-White', (0., 0., 1.)), ('Green-White-linear', (0., 1., 0.)), ('Green-White-exponential', (0., 1., 0.)), ] '''predefined list of colormaps suitable for anatomical volumes''' anatomicalFusionColormaps = [('B-W LINEAR-fusion', (1., 1., 1.)), ('Blue-White-fusion', (0., 0., 1.)), ('Green-White-linear-fusion', (0., 1., 0.)), ] '''predefined list of colormaps suitable for fusionned anatomical volumes''' diffusionColormaps = anatomicalColormaps '''predefined list of colormaps suitable for diffusion volumes''' diffusionFusionColormaps = anatomicalFusionColormaps '''predefined list of colormaps suitable for fusionned diffusion volumes''' functionalColormaps = [('RED TEMPERATURE', (1., 0.5, 0.)), ('RAINBOW', (1., 0., 0.)), ('Blue-Red', (1., 0., 0.)), ('actif-ret', (1., 1., 0.)), ('Yellow-red', (1., 1., 0.)), ] '''predefined list of colormaps suitable for functional volumes''' functionalFusionColormaps = [('Rainbow1-fusion', (1., 0., 0.)), ('Blue-Red-fusion', (1., 0., 0.)), ('Yellow-red-fusion', (1., 1., 0.)), ] '''predefined list of colormaps suitable for fusionned functional volumes''' twotailColormaps = [('tvalues100-200-100-lfusion', (1., 0., 0.)), ('tvalues100-100-100-lfusion', (1., 0., 0.)), ] '''predefined list of colormaps suitable for two-tail T-values volumes''' twotailFusionColormaps = [('tvalues100-200-100', (1., 0., 0.)), ('tvalues100-100-100', (1., 0., 0.)), ] '''predefined list of colormaps suitable for fusionned two-tail T-values volumes''' labelsColormaps = [('Blue-Red', (1., 0., 0.)), ('Talairach', (0., .0, .0, )), ] '''predefined list of colormaps suitable for labels volumes''' labelsFusionColormaps = [ ] '''predefined list of colormaps suitable for fusionned labels volumes''' binaryColormaps = [('BLUE-lfusion', (0., 0., 1.)), ('GREEN-lfusion', (0., 1., 0.)), ('RED-lfusion', (1., 0., 0.)), ('CYAN-lfusion', (0., 1., 1.)), ('VIOLET-lfusion', (1., 0., 1.)), ('YELLOW-lfusion', (1., 1., 0.)), ('WHITE-lfusion', (1., 1., 1.)), ] '''predefined list of colormaps suitable for binary volumes''' binaryFusionColormaps = [('BLUE-ufusion', (0., 0., 1.)), ('GREEN-ufusion', (0., 1., 0.)), ('RED-ufusion', (1., 0., 0.)), ('CYAN-ufusion', (0., 1., 1.)), ('VIOLET-ufusion', (1., 0., 1.)), ('YELLOW-ufusion', (1., 1., 0.)), ('Black-ufusion', (1., 1., 1.)), ] '''predefined list of colormaps suitable for fusionned binary volumes'''
[docs]def chooseColormaps(vols): '''Automatically chooses distinc colormaps for a list of volumes - returns: a list of colormaps names. They should be known from Anatomist. ''' cmapbytype = [anatomicalColormaps, diffusionColormaps, functionalColormaps, labelsColormaps, twotailColormaps, binaryColormaps] cmapbytype_fusion = [anatomicalFusionColormaps, diffusionFusionColormaps, functionalFusionColormaps, labelsFusionColormaps, twotailColormaps, binaryColormaps] def findRgbColormap(cmap): for cmt in cmapbytype: for cm in cmt: if cm[0] == cmap: return cm[1] return None chooseColormaps.findRgbColormap = findRgbColormap def rgbDist(rgb1, rgb2): rgb = numpy.array(rgb1) - rgb2 return numpy.sqrt((rgb ** 2).sum()) chooseColormaps.rgbDist = rgbDist def chooseColormap(vtype, cmaps, rgbcmaps): disttoothers = [] ct = cmapbytype[vtype] for cmap in ct: if cmap[0] in cmaps: disttoothers.append(0.) continue # find min distance from cmap to already used ones bd = numpy.sqrt(3.) for c in rgbcmaps: if c is not None: dist = rgbDist(cmap[1], c) if dist < bd: bd = dist disttoothers.append(bd) disttoothers = numpy.array(disttoothers) bests = numpy.where(disttoothers == disttoothers.max())[0] if len(bests) == 1: return ct[bests[0]] else: candidates = numpy.array([ct[x][0] for x in bests]) count = numpy.zeros(len(bests)) for c in cmaps: cc = numpy.where(candidates == c)[0] if len(cc) == 1: count[cc[0]] += 1 return ct[bests[count.argmin()]] cmaps = [] rgbcmaps = [] typeclasses = [{}, {}, {}, {}, {}, {}] vtypes = [] i = 0 for vol in vols: lh = 1. if type(vol) is not type({}): vol = checkVolume(vol) if vol['binary_volume']: vtype = 5 else: vtype = numpy.argmax(vol['volume_contents_likelihoods']) if vtype == 2 and vol['negative_values']: vtype = 4 else: # classify by likelihood lh = vol['volume_contents_likelihoods'][vtype] tc = typeclasses[vtype] if lh in tc: tc[lh].append((vol, i)) else: tc[lh] = [(vol, i)] vtypes.append(vtype) i += 1 # re-order volumes by type/likelihood ordered = [] for t in typeclasses: keys = list(t.keys()) keys.sort() keys.reverse() # highest at first for k in keys: ordered += t[k] # now choose cmap in ordered volumes for vol, i in ordered: pal = vol.get('palette') if pal: if pal not in cmaps: cmaps.append(pal) rgbcmaps.append(findRgbColormap(pal)) else: pal = None if not pal: vtype = vtypes[i] cmap = chooseColormap(vtype, cmaps, rgbcmaps) cmaps.append(cmap[0]) rgbcmaps.append(cmap[1]) orderedcmaps = [None] * len(ordered) for i in range(len(ordered)): orderedcmaps[ordered[i][1]] = cmaps[i] return orderedcmaps