Source code for ClearMap.ImageProcessing.Topology.Topology2d

# -*- coding: utf-8 -*-
"""
Topology2d
==========

Defines basic 2d discrete topology utils.
"""
__author__    = 'Christoph Kirst <christoph.kirst.ck@gmail.com>'
__license__   = 'GPLv3 - GNU General Pulic License v3 (see LICENSE.txt)'
__copyright__ = 'Copyright © 2020 by Christoph Kirst'
__webpage__   = 'http://idisco.info'
__download__  = 'http://www.github.com/ChristophKirst/ClearMap2'



import numpy as np

###############################################################################
### Topology
###############################################################################

# Topology numbers T8 and T4bar in 2D
#Configuration:
# 3 2 1			
# 4 X 0
# 5 6 7

t4_bar = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1,
       1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 2, 2,
       2, 2, 1, 2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 1, 2,
       1, 2, 1, 2, 2, 3, 2, 2, 1, 2, 1, 1, 1, 2, 1, 1, 0, 1, 0, 1, 1, 2, 1,
       1, 0, 1, 0, 1, 1, 2, 1, 2, 1, 2, 1, 2, 2, 3, 2, 2, 1, 2, 1, 1, 1, 2,
       1, 1, 0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0, 1, 1, 2, 1, 2, 1, 2, 1, 2, 2,
       3, 2, 2, 1, 2, 1, 2, 2, 3, 2, 2, 1, 2, 1, 2, 2, 3, 2, 2, 1, 2, 1, 2,
       2, 3, 2, 3, 2, 3, 2, 3, 3, 4, 3, 3, 2, 3, 2, 2, 2, 3, 2, 2, 1, 2, 1,
       2, 2, 3, 2, 2, 1, 2, 1, 1, 1, 2, 1, 2, 1, 2, 1, 2, 2, 3, 2, 2, 1, 2,
       1, 1, 1, 2, 1, 1, 0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0, 1, 1, 2, 1, 2, 1,
       2, 1, 2, 2, 3, 2, 2, 1, 2, 1, 1, 1, 2, 1, 1, 0, 1, 0, 1, 1, 2, 1, 1,
       0, 1, 0]);
  
t8 = np.array([0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1,
       1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 2, 2,
       2, 2, 1, 2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 1, 2,
       1, 2, 1, 2, 2, 3, 2, 2, 1, 2, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 2, 1,
       1, 1, 1, 1, 1, 1, 2, 1, 2, 1, 2, 1, 2, 2, 3, 2, 2, 1, 2, 1, 1, 1, 2,
       1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 2, 1, 2, 1, 2, 1, 2, 2,
       3, 2, 2, 1, 2, 1, 2, 2, 3, 2, 2, 1, 2, 1, 2, 2, 3, 2, 2, 1, 2, 1, 2,
       2, 3, 2, 3, 2, 3, 2, 3, 3, 4, 3, 3, 2, 3, 2, 2, 2, 3, 2, 2, 1, 2, 1,
       2, 2, 3, 2, 2, 1, 2, 1, 1, 1, 2, 1, 2, 1, 2, 1, 2, 2, 3, 2, 2, 1, 2,
       1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 2, 1, 2, 1,
       2, 1, 2, 2, 3, 2, 2, 1, 2, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1,
       1, 1, 1])


###############################################################################
### Labels and Indices
###############################################################################


[docs]def plane_label(center = None):
  if center is not None:
    base = np.array([[3,2,1],[4,0,0],[5,6,7]]);
    base = np.power(2, base);
    base[1,1] = center;
  else:
    base = np.array([[4,3,2],[5,0,1],[6,7,8]]);
    base = np.power(2, base);
  return base;




[docs]def plane_base_2(center = None):
  """Returns an array with base 2 numbers on the plane for convolution and lut matching"""
  plane = np.zeros((3,3), dtype = int);
  k = 0;
  for y in range(3):
    for x in range(3):
      if center is not None and x == 1 and y ==1:
        plane[x,y] = center;
      else:
        plane[x,y] = 2**k;
        k+=1;
  return plane;




[docs]def plane_from_index(index, center = None):
  """Returns a boolean plane for the corresponding index"""
  plane = np.zeros((3,3), dtype = bool);
  d = 0;
  for y in range(3):
    for x in range(3):
      if center is not None and x == 1 and y == 1:
        plane[x,y] = center;
      else:
        plane[x,y] = (index >> d) & 0x01;
        d += 1;
  return plane;


  

[docs]def plane_to_index(plane, center = None):
  """Returns index for a boolean cube"""
  return (plane_base_2(center=center) * np.array(plane)).sum()



###############################################################################
### Neighbourhoods
###############################################################################

#TODO: still needed ? -> clean up

[docs]def extract_neighbourhood(img,x,y):
  """Return the neighbourhoods of the indicated voxels
  
  Arguments:
    img (array): the 2d image
    x,y (n array): coordinates of the voxels to extract neighbourhoods from
  
  Returns:
    array (nx9 array): neighbourhoods
    
  Note:
    Assumes borders of the image are zero so that 0<x,y<w,h !
  """
  nhood = np.zeros((x.shape[0],9), dtype = bool);
  
  # calculate indices (if many voxels this is only 9 loops!)
  for xx in range(3):
    for yy in range(3):
        #w = _xyz_to_neighbourhood[xx,yy,zz];
        w = 3 * xx + yy;
        idx = x+xx-1; idy = y+yy-1;
        nhood[:,w]=img[idx, idy];
  
  nhood.shape = (nhood.shape[0], 3, 3);
  nhood[:, 1, 1] = 0;
  return nhood;



###############################################################################
### Testing
###############################################################################

def _test():
  import ClearMap.ImageProcessing.Topology.Topology2d as t2d;
  
  p = [[1,0,0],[0,0,1],[1,0,0]];
  pi = t2d.plane_to_index(p)
  print('t4b: %d, t8: %d' % (t2d.t4_bar[pi], t2d.t8[pi]))
  
  p = [[1,1,0],[0,0,1],[1,0,0]];
  pi = t2d.plane_to_index(p)
  print('t4b: %d, t8: %d' % (t2d.t4_bar[pi], t2d.t8[pi]))

  p2 = t2d.plane_from_index(pi)
  assert np.all(p == p2)