Source code for ClearMap.ImageProcessing.Tracing.Connect

# -*- coding: utf-8 -*-
"""
Module to Postprocess skeletonized vasculature data
"""
__author__    = 'Christoph Kirst <ckirst@rockefeller.edu>'
__license__   = 'MIT License <http://www.opensource.org/licenses/mit-license.php>'
__copyright__ = 'Copyright 2017 by Christoph Kirst, The Rockefeller University, New York City'


import numpy as np

import scipy.ndimage as ndi

#TODO:
import ClearMap.ParallelProcessing.DataProcessing.ArrayProcessing as ap
import ClearMap.ParallelProcessing.DataProcessing.ConvolvePointList as cpl

import ClearMap.ImageProcessing.Tracing.Trace as trc
import ClearMap.ImageProcessing.Differentiation.Hessian as cur

import ClearMap.ImageProcessing.Topology.Topology3d as t3d

import ClearMap.Visualization.Plot3d as p3d



[docs]
def extractNeighbourhood(data, center, radius):
  """Extract local neighborhood with specific radius, if to close to border pad with zeros"""
  order = 'C';
  if np.isfortran(data):
    order = 'F';
  
  if isinstance(radius, int):
    radius = (radius,);
  radius = radius * 3;
  radius = np.array(radius[:3]);
  
  dlo = [max(0, c-r) for c,r in zip(center, radius)];
  dhi = [min(s, c+r+1) for c,r,s in zip(center, radius, data.shape)];
  
  nlo = [-min(0, c-r) for c,r in zip(center, radius)];
  nhi = [2*r+1 + min(0, s-c-r-1) for c,r,s in zip(center, radius, data.shape)];
  
  #print center
  #print radius
  #print dlo, dhi
  #print nlo, nhi  
  
  nbh = np.zeros(2*radius + 1, dtype = data.dtype, order = order);
  nbh[nlo[0]:nhi[0],nlo[1]:nhi[1],nlo[2]:nhi[2]] = data[dlo[0]:dhi[0],dlo[1]:dhi[1],dlo[2]:dhi[2]];
  
  return nbh;




[docs]
def findEndpoints(skel, points, border = None):
  """Find endpoints in skeleton to try to reconnect"""
  order = 'C';
  if np.isfortran(skel):
    order = 'F'; 
  
  #find node degrees
  deg = cpl.convolve_3d_indices(skel, t3d.n26, points, out_dtype = 'uint8')
  
  #isolated and end points
  isolated = points[deg == 0];
  ends = points[deg == 1];
  
  # remove border points
  if border is not None:
    if not isinstance(border, tuple):
      border = (border,);
    border = border * 3;
    border = border[:3];
  
    x,y,z = np.unravel_index(isolated, skel.shape, order = order)
    ids = np.ones(isolated.shape, dtype = bool);
    for xx,ss,bb in zip((x,y,z), skel.shape, border):
      ids = np.logical_and(ids, np.logical_and(bb < xx, xx < ss - bb));
    isolated = isolated[ids];
    
    x,y,z = np.unravel_index(ends, skel.shape, order = order)
    ids = np.ones(ends.shape, dtype = bool);
    for xx,ss,bb in zip((x,y,z), skel.shape, border):
      ids = np.logical_and(ids, np.logical_and(bb < xx, xx < ss - bb));
    ends = ends[ids];
  
  return ends, isolated




[docs]
def addPathToMask(mask, path, value = True):
  order = 'C';
  if np.isfortran(mask):
    order = 'F'; 
  m = np.reshape(mask, -1, order = order);
  ids = np.ravel_multi_index(path.T, mask.shape, order = order);
  m[ids] = value;

  #return mask;



[docs]
def addDilatedPathToMask(mask, path, iterations = 1):
  order = 'C';
  if np.isfortran(mask):
    order = 'F'; 
  ids = np.ravel_multi_index(path.T, mask.shape, order = order);
  m = np.zeros(mask.shape, dtype = bool, order = order);
  mf = np.reshape(m, -1, order = order);
  mf[ids] = True;
  m = ndi.morphology.binary_dilation(m, structure = np.ones((3,3,3), dtype = bool), iterations = iterations);
  m = np.asarray(m, order = order);

  #return np.logical_or(mask, m);


viewer = None;


[docs]
def plotData(data, skel, binary, ends = None, isolated = None, replot = True):
  global viewer
  
  img = np.zeros(data.shape, dtype = int);
  img[:] = img + binary + skel;
  
  if ends is not None:
    xyz = np.array(np.unravel_index(ends, skel.shape, order = 'C')).T;
    for x,y,z in xyz:
      img[x,y,z] = 3;
  
  if isolated is not None:  
    xyz = np.array(np.unravel_index(isolated, skel.shape, order = 'C')).T;
    for x,y,z in xyz:
      img[x,y,z] = 4;
  
  if replot:
    try:
      viewer[0].setSource(data);
      viewer[1].setSource(data);
    except:
      viewer = p3d.plot([data, img]);
  else:
    return p3d.plot([data, img]);


 

[docs]
def plotTracingResult(path, data_nbh, mask_nbh, center, radius, tubeness_nbh, skeleton = None, distance_nbh = None):
  global viewer;
  
  img_nbh = np.asarray(data_nbh, dtype = float);
  bin_nbh = np.zeros(mask_nbh.shape, dtype = int);   
  
  if skeleton is not None:
    bin_nbh += extractNeighbourhood(skeleton, center, radius); 

  addPathToMask(bin_nbh, path, value = 4);
  if points is not None:
    pts = points - center + radius;
    for d in range(3):
      ids = np.logical_and(0 <= pts[:,d], pts[:,d] < data_nbh.shape[d]);
      pts = pts[ids];
    addPathToMask(bin_nbh, pts, value = 3);
  
  addPathToMask(img_nbh, path, value = 512);
  addPathToMask(tubeness_nbh, path, value = 50);
  if distance_nbh is not None:
    addPathToMask(distance_nbh, path, value = 50);
  
  try:
    viewer[0].setSource(bin_nbh);
    viewer[1].setSource(img_nbh);
    viewer[2].setSource(tubeness_nbh);
    if distance_nbh is not None:
      viewer[3].setSource(distance_nbh);
  except:
    if distance_nbh is not None:
      viewer = p3d.plot([bin_nbh, img_nbh, tubeness_nbh, distance_nbh])
    else:
      viewer = p3d.plot([bin_nbh, img_nbh, tubeness_nbh])










[docs]
def tracePointToMask(data, mask, center, radius, points = None, plot = False, skeleton = None, tubeness = None, 
                     removeLocalMask = True,
                     maxSteps = 500, verbose = False, **trace_parameter):
  """Trace an endpoint to a mask"""
    
  # cut out neighbourhood from data
  center_nbh = np.array(center) - center + radius; 
  data_nbh = extractNeighbourhood(data, center, radius);
  mask_nbh = extractNeighbourhood(mask, center, radius);
  
  if tubeness is None:
    tubeness_nbh = cur.tubeness(ndi.gaussian_filter(np.asarray(data_nbh, dtype = float), sigma = 1.0));
  else:
    tubeness_nbh = extractNeighbourhood(tubeness, center, radius);
  
  
  if removeLocalMask is not None:
    mask_nbh_label, mask_label = ndi.label(mask_nbh, structure = np.ones((3,3,3), dtype = bool));
    ids = mask_nbh_label[tuple(center_nbh)] == mask_nbh_label;
    mask_nbh[ids] = False;
  
  distance_nbh = ndi.distance_transform_edt(np.logical_not(mask_nbh))
  
  path, quality = trc.traceToMask(np.asarray(data_nbh, dtype = float), tubeness_nbh, center_nbh, distance_nbh, 
                                  maxSteps = maxSteps, verbose = verbose, returnQuality = True, **trace_parameter);
  
  if verbose:
    if len(path) > 0:
      print('Path of length = %d with quality = %f (per length = %f)' % (len(path), quality, quality / len(path)));
    else:
      print('No path found!');
  
  if plot:
    plotTracingResult(path, data_nbh, mask_nbh, center, radius, tubeness_nbh, skeleton = skeleton, distance_nbh = distance_nbh);
  
  return path + center - center_nbh, quality;




[docs]
def tracePointToNeighbor(data, mask, center, neighbor, radius, points = None, plot = False, skeleton = None, tubeness = None, 
                         removeLocalMask = True,
                         maxSteps = 500, verbose = False, **trace_parameter):
  """Trace an endpoint to a neighbour"""
  # cut out neighbourhood from data
  center_nbh = np.array(center) - center + radius; 
  neighbor_nbh = np.array(neighbor) - center + radius;   
  
  data_nbh = extractNeighbourhood(data, center, radius);
  mask_nbh = extractNeighbourhood(mask, center, radius);
  
  if tubeness is None:
    tubeness_nbh = cur.tubeness(ndi.gaussian_filter(np.asarray(data_nbh, dtype = float), sigma = 1.0));
  else:
    tubeness_nbh = extractNeighbourhood(tubeness, center, radius);
  
  if removeLocalMask is not None:
    mask_nbh_label, mask_label = ndi.label(mask_nbh, structure = np.ones((3,3,3), dtype = bool));
    ids = mask_nbh_label[tuple(center_nbh)] == mask_nbh_label;
    mask_nbh[ids] = False;
  
  #distance_nbh = ndi.distance_transform_edt(np.logical_not(mask_nbh))
  
  path, quality = trc.trace(np.asarray(data_nbh, dtype = float), tubeness_nbh, 
                   center_nbh, neighbor_nbh, 
                   maxSteps = maxSteps, verbose = verbose, returnQuality = True, **trace_parameter);
  
  if verbose:
    if len(path) > 0:
      print('Path of length = %d with quality = %f (per length = %f)' % (len(path), quality, quality / len(path)));
    else:
      print('No path found!');
  
  if plot:
    plotTracingResult(path, data_nbh, mask_nbh, center, radius, tubeness_nbh, skeleton = skeleton);
  
  return path + center - center_nbh, quality;




[docs]
def connectPoint(data, mask, endpoints, start_index, radius, 
                 tubeness = None, min_quality = None, remove_local_mask = True,
                 skeleton = None,
                 verbose = False, **trace_parameter):
  """Tries to connect an end point"""
  
  #outine:
  # find neighbour end points and try to connect to nearest one
  # if path score good enough add path and remove two endpoints
  # else try to connect to binarized image
  # if path score good enugh connect to closest skeleton point
  # else not connectable
  
  #assumes everything is in fotran order
  strides = np.array(data.strides) / data.itemsize;
  shape = data.shape;
  #print strides, shape
  
  center_flat = endpoints[start_index];
  center_xyz  = np.array(np.unravel_index(center_flat, data.shape, order = 'F'));
  
  mask_nbh = extractNeighbourhood(mask, center_xyz, radius);
  data_nbh = np.asarray(extractNeighbourhood(data, center_xyz, radius), dtype = float, order = 'F');
  shape_nbh = mask_nbh.shape;
  
  center_nbh_xyz = np.zeros(3, dtype = int) + radius; 
  #center_nbh_flat = np.ravel_multi_index(center_nbh_xyz, shape_nbh, order = 'F');
  
  if tubeness is None:
    tubeness_nbh = cur.tubeness(ndi.gaussian_filter(np.asarray(data_nbh, dtype = float), sigma = 1.0));
    tubeness_nbh = np.asarray(tubeness_nbh, order = 'F');
  else:
    tubeness_nbh = extractNeighbourhood(tubeness, center_xyz, radius);
  
  mask_nbh_label = np.empty(shape_nbh, dtype = 'int32', order = 'F');
  _ = ndi.label(mask_nbh, structure = np.ones((3,3,3), dtype = bool), output = mask_nbh_label);
  local_nbh = mask_nbh_label[tuple(center_nbh_xyz)] == mask_nbh_label;  
  
  # end point neighbours
  nbs_flat = ap.findNeighbours(endpoints, start_index, shape, strides, radius);
  
  if len(nbs_flat) > 0:
    nbs_nbh_xyz  = np.vstack(np.unravel_index(nbs_flat, shape, order = 'F')).T - center_xyz + center_nbh_xyz;
    nbs_nbh_flat = np.ravel_multi_index(nbs_nbh_xyz.T, shape_nbh, order = 'F');
    
    # remove connected neighbours
    non_local_nbh_flat = np.reshape(np.logical_not(local_nbh), -1, order = 'F');
    nbs_nbh_non_local_flat = nbs_nbh_flat[non_local_nbh_flat[nbs_nbh_flat]];
    
    if len(nbs_nbh_non_local_flat) > 0:
      #find nearest neighbour
      nbs_nbh_non_local_xyz  = np.vstack(np.unravel_index(nbs_nbh_non_local_flat, shape, order = 'F')).T;
      
      nbs_nbh_non_local_dist = nbs_nbh_non_local_xyz - center_nbh_xyz;
      nbs_nbh_non_local_dist = np.sum(nbs_nbh_non_local_dist * nbs_nbh_non_local_dist, axis = 1);
      
      neighbor_nbh_xyz = nbs_nbh_non_local_xyz[np.argmin(nbs_nbh_non_local_dist)];
      
      path, quality = trc.trace(data_nbh, tubeness_nbh, 
                                center_nbh_xyz, neighbor_nbh_xyz, 
                                verbose = False, returnQuality = True, **trace_parameter);
      
      if len(path) > 0:
        if quality / len(path) < min_quality:
          if verbose:
            print('Found good path to neighbour of length = %d with quality = %f (per length = %f) [%d / %d nonlocal neighbours]' % (len(path), quality, quality / len(path), len(nbs_nbh_non_local_flat), len(nbs_flat)));  #print path
          return path + center_xyz - center_nbh_xyz, quality;
        else:
          if verbose:
            print('Found bad  path to neighbour of length = %d with quality = %f (per length = %f) [%d / %d nonlocal neighbours]' % (len(path), quality, quality / len(path), len(nbs_nbh_non_local_flat), len(nbs_flat)));  #print path
      else:
        if verbose:
          print('Found no path to neighbour [%d / %d nonlocal neighbours]' % (len(nbs_nbh_non_local_flat), len(nbs_flat)));  #print path
    
    
  #tracing to neares neighbour failed    
  if verbose:
    print('Found no valid path to neighbour, now tracing to binary!');  #print path

  # Tracing to next binary
  if remove_local_mask:
    mask_nbh[local_nbh] = False; 
  
  distance_nbh = ndi.distance_transform_edt(np.logical_not(mask_nbh))
  distance_nbh = np.asarray(distance_nbh, order = 'F');
  
  path, quality = trc.traceToMask(data_nbh, tubeness_nbh, center_nbh_xyz, distance_nbh, 
                                  verbose = False, returnQuality = True, **trace_parameter);
  
  if len(path) > 0:
    if quality / len(path) < min_quality:
      if verbose:
        print('Found good path to binary of length = %d with quality = %f (per length = %f)' % (len(path), quality, quality / len(path)));  #print path
            
      # trace to skeleton
      if skeleton is not None:
        #find closest point on skeleton
        final_xyz = path[0];
        skeleton_nbh = extractNeighbourhood(skeleton, center_xyz, radius);
        local_end_path_nbh = mask_nbh_label[tuple(final_xyz)] == mask_nbh_label;
        skeleton_nbh_dxyz = np.vstack(np.where(np.logical_and(skeleton_nbh, local_end_path_nbh))).T - final_xyz;
        if len(skeleton_nbh_dxyz) == 0: # could not find skeleton nearby -> give up for now
          return path + center_xyz - center_nbh_xyz, quality;
        
        skeleton_nbh_dist = np.sum(skeleton_nbh_dxyz * skeleton_nbh_dxyz, axis = 1);
        closest_dxyz =  skeleton_nbh_dxyz[np.argmin(skeleton_nbh_dist)];
        closest_xyz = closest_dxyz + final_xyz;
        #print path[0], path[-1]
        #print center_nbh_xyz, closest_dxyz
        
        #generate pixel path
        max_l = np.max(np.abs(closest_dxyz)) + 1;
        path_add_xyz  = np.vstack([np.asarray(np.linspace(f, c, max_l), dtype = int) for f,c in zip(final_xyz, closest_xyz)]).T;
        path_add_flat = np.ravel_multi_index(path_add_xyz.T, shape_nbh);
        _, ids = np.unique(path_add_flat, return_index = True);
        path_add_xyz = path_add_xyz[ids];
        #print path_add_xyz;
        path = np.vstack([path, path_add_xyz]); # note: this is not an ordered path anymore!
      
      return path + center_xyz - center_nbh_xyz, quality;
    else:
      if verbose:
        print('Found bad  path to binary of length = %d with quality = %f (per length = %f)' % (len(path), quality, quality / len(path)));  #print path
   
  if verbose:
    print('Found no valid path to binary!');
    
  return np.zeros((0,3)), 0;



import ClearMap.ParallelProcessing.SharedMemoryManager as smm;
#import ClearMap.ParallelProcessing.SharedMemoryProcessing as smp;

import multiprocessing as mp

import tempfile as tmpf



[docs]
def order(array):
  if np.isfortran(array):
    return 'F';
  else:
    return 'C';




[docs]
def processSingleConnection(args):
  global temporary_folder;
  i = args;
  
  #steps_done += 1;
  #if steps_done % 100 == 0:
  #  print('Processing step %d / %d...' % (steps_done, steps_total));
  if i % 100 == 0:
    print('Processing step %d...' % i);

  try:  
    fid = mp.current_process()._identity[0];
  except: # mostlikely sequential mode
    fid = 0;
  
  data = smm.get(0);
  mask = smm.get(1);
  skel = smm.get(2);
  spts = smm.get(3);
  
  res = connectPoint(data, mask, spts, i, skeleton = skel, 
                     radius = 20, 
                     tubeness = None, remove_local_mask = True, 
                     min_quality = 15.0,
                     verbose = False,
                     maxSteps = 12000, costPerDistance = 1.0);
  #return res;
  
  path, score = res;
  if len(path) > 0: 
    #convert path to indices
    path_flat = np.ravel_multi_index(path.T, data.shape, order = order(data));
     
    #append to file
    try:
      fid = mp.current_process()._identity[0];
    except:
      fid = 0; # sequential mode
    fn  = '%s/path_%03d.tmp' % (temporary_folder, fid);
    #print fn;
    
    fb = open(fn,"ab");
    path_flat.tofile(fb);
    fb.close();                     
   
  return;



import os;
import shutil;
import ClearMap.Utils.Timer as tmr

import gc

temporary_folder = None;  


[docs]
def addConnections(data, mask, skeleton, points, radius = 20, 
                   start_points = None, 
                   remove_local_mask = True, min_quality = 15.0,
                   add_to_skeleton = True, add_to_mask = False,
                   verbose = True, processes = mp.cpu_count(), block_size = 5000, debug = False):
  global temporary_folder;
  
  timer = tmr.Timer();
  timer_total = tmr.Timer();
  
  if start_points is None:
    ends, isolated = findEndpoints(skeleton, points, border = 20);
    start_points = np.hstack([ends, isolated]);
    start_points = smm.asShared(start_points);
  npts = len(start_points);
  
  if verbose:
    timer.printElapsedTime('Found %d endpoints' % (npts,));
    timer.reset();
  
  assert smm.isShared(data);
  assert smm.isShared(mask);
  assert smm.isShared(skeleton);
  #steps_total = npts / processes;
  #steps_done = 0;
  
  smm.clean();
  data_hdl = smm.insert(data);
  mask_hdl = smm.insert(mask);
  skel_hdl = smm.insert(skeleton);
  spts_hdl = smm.insert(start_points);
  if not (data_hdl == 0 and mask_hdl == 1 and skel_hdl == 2 and spts_hdl == 3):
    raise RuntimeError('The shared memory handles are invalid: %d, %d, %d, %d' % (data_hdl, mask_hdl, skel_hdl, spts_hdl));
  
  #generate temporary folder to write path too
  temporary_folder = tmpf.mkdtemp();
  
  if verbose:
    timer.printElapsedTime('Preparation of %d connections' % (npts,));
    timer.reset();
    
  # process in parallel / block processing is to clean up memory leaks for now
  nblocks = max(1, int(np.ceil(1.0 * npts / processes / block_size)));
  ranges = np.asarray(np.linspace(0, npts, nblocks+1), dtype = int);
  #nblocks = 1;
  #ranges = [0, 100];
  for b in range(nblocks):
    argdata = np.arange(ranges[b], ranges[b+1]);
    if debug:
      result = [processSingleConnection(a) for a in argdata];
    else:
      #mp.process.current_process()._counter = mp.process.itertools.count(1); #reset worker counter
      pool = mp.Pool(processes = processes);
      #for i,_ in enumerate(pool.imap_unordered(processSingleConnection, argdata)):
      #  if i % 100 == 0:
      #    timer.printElapsedTime('Iteration %d / %d' % (i + ranges[b], npts));
      pool.map(processSingleConnection, argdata)
          
      pool.close();
      pool.join();
    gc.collect();
  
  smm.free(data_hdl);
  smm.free(mask_hdl);
  smm.free(skel_hdl);
  smm.free(spts_hdl);
  
  if verbose:
    timer.printElapsedTime('Processing of %d connections' % (npts,));
    timer.reset();

  #return result;

  #get paths from temporay files
  result = [];
  for f in os.listdir(temporary_folder):
    result.append(np.fromfile(os.path.join(temporary_folder,f), dtype = int));
  result = np.hstack(result);
  
  #clean up
  shutil.rmtree(temporary_folder);
  temporary_folder = None;
  
  if verbose:
    timer.printElapsedTime('Loading paths with %d points' % (len(result),));
    timer.reset();
  
  #add dilated version to skeleton
  if add_to_skeleton:
    skeleton_f = np.reshape(skeleton, -1, order = 'A');
    strides = skeleton.strides;
    skeleton_f_len = len(skeleton_f);
    for dx in [-1,0,1]:
      for dy in [-1,0,1]:
        for dz in [-1,0,1]:
          offset= dx * strides[0] + dy * strides[1] + dz * strides[2];
          pos = result + offset;
          pos = pos[np.logical_and(pos >= 0, pos < skeleton_f_len)];
          skeleton_f[pos] = True;
  
    if verbose:
      timer.printElapsedTime('Added paths with %d points to skeleton' % (len(result),));
      #add dilated version to skeleton
  
  if add_to_mask:
    mask_f = np.reshape(mask, -1, order = 'A');
    strides = mask.strides;
    for dx in [-1,0,1]:
      for dy in [-1,0,1]:
        for dz in [-1,0,1]:
          offset= dx * strides[0] + dy * strides[1] + dz * strides[2];
          pos = result + offset;
          pos = pos[np.logical_and(pos >= 0, pos < skeleton_f_len)];
          mask_f[pos] = True;
  
    if verbose:
      timer.printElapsedTime('Added paths with %d points to mask' % (len(result),));  
      
      
    timer_total.printElapsedTime('Connection post-processing added %d points' % (len(result),));
  
  return result;



def _test():
  #%%
  import numpy as np
  import scipy.ndimage as ndi
  import ClearMap.DataProcessing.LargeData as ld
  import ClearMap.Visualization.Plot3d as p3d
  import ClearMap.DataProcessing.ConvolvePointList as cpl
  import ClearMap.ImageProcessing.Skeletonization.Topology3d as t3d
  import ClearMap.ImageProcessing.Skeletonization.SkeletonCleanUp as scu
  
  import ClearMap.ImageProcessing.Tracing.Connect as con;
  reload(con);  
  
  data   = np.load('/home/ckirst/Desktop/data.npy');
  binary = np.load('/home/ckirst/Desktop/binarized.npy');
  skel   = np.load('/home/ckirst/Desktop/skel.npy');
  #points = np.load('/home/ckirst/Desktop/pts.npy');
  
  data   = np.copy(data,    order = 'F');
  binary = np.copy(binary,  order = 'F');
  skel   = np.copy(skel,    order = 'F');
  skel_copy = np.copy(skel, order = 'F');
  points = np.ravel_multi_index(np.where(skel), skel.shape, order = 'F');
  
  skel, points = scu.cleanOpenBranches(skel, skel_copy, points, length = 3, clean = True);
  deg = cpl.convolve3DIndex(skel, t3d.n26, points);
  
  ends, isolated = con.findEndpoints(skel, points, border = 25);
  special = np.sort(np.hstack([ends, isolated]))
  
  ends_xyz = np.array(np.unravel_index(ends, data.shape, order = 'F')).T;
  isolated_xyz = np.array(np.unravel_index(isolated, data.shape, order = 'F')).T;
  special_xyz = np.vstack([ends_xyz, isolated_xyz]);
  
  
  #%%
  import ClearMap.ParallelProcessing.SharedMemoryManager as smm
  data_s   = smm.asShared(data, order = 'F');
  binary_s = smm.asShared(binary.view('uint8'), order = 'F');
  skel_s   = smm.asShared(skel.view('uint8'),   order = 'F');
  
  smm.clean();
  res = con.addConnections(data_s, binary_s, skel_s, points, radius = 20, 
                           start_points = None,
                           add_to_skeleton=True, add_to_mask = True,
                           verbose = True, processes = 4, debug = False, block_size=10)
  
  skel_s = skel_s.view(bool);
  binary_s = binary_s.view(bool);

  #%%
  mask_img = np.asarray(binary, dtype= int, order = 'A'); 
  mask_img[:] = mask_img + binary_s;
  mask_img[:] = mask_img + skel;
  
  data_img = np.copy(data, order = 'A');
  data_img[skel] = 120;
  
  mask_img_f = np.reshape(mask_img, -1, order = 'A');
  data_img_f = np.reshape(data_img, -1, order = 'A');
  
  mask_img_f[res] = 7;
  data_img_f[res] = 512;
  
  mask_img_f[special] = 8;
  data_img_f[special] = 150;
  
  for d in [3,4,5]:
    mask_img_f[points[deg == d]] = d+1;
  
  try:
    con.viewer[0].setSource(mask_img);
    con.viewer[1].setSource(data_img);
  except:
    con.viewer = p3d.plot([mask_img, data_img]);
  
  con.viewer[0].setMinMax([0, 8]);
  con.viewer[1].setMinMax([24, 160]);          

  
  #%%
  mask = binary; 
  data_new = np.copy(data, order = 'A');
  data_new[skel] = 120;
  
  skel_new = np.asarray(skel, dtype= int, order = 'A');
  skel_new[:] = skel_new + binary;
  
  binary_new = np.copy(binary, order = 'A');
  qs = [];
  for i,e in enumerate(special):
    print('------');
    print('%d / %d'  % (i, len(special)));
    path, quality = con.connectPoint(data, mask, special, i, radius = 25,
                                     skeleton = skel, 
                                     tubeness = None, remove_local_mask = True, 
                                     min_quality = 15.0,
                                     verbose = True,
                                     maxSteps = 15000, costPerDistance = 1.0);

    
    #print path, quality
    if len(path) > 0:
      qs.append(quality * 1.0/len(path));      
      
      q = con.addPathToMask(skel_new, path, value = 7);
      q = con.addPathToMask(data_new, path, value = 512);
      binary_new = con.addDilatedPathToMask(binary_new, path, iterations = 1);
 
  skel_new[:] = skel_new + binary_new;
  q = con.addPathToMask(skel_new, special_xyz, value = 6);
  for d in [3,4,5]:
    xyz = np.array(np.unravel_index(points[deg == d], data.shape, order = 'F')).T;
    q = con.addPathToMask(skel_new, xyz, value = d);  
  q = con.addPathToMask(data_new, special_xyz, value = 150);
    
  try:
    con.viewer[0].setSource(skel_new);
    con.viewer[1].setSource(data_new);
  except:
    con.viewer = p3d.plot([skel_new, data_new]);
  
  con.viewer[0].setMinMax([0, 8]);
  con.viewer[1].setMinMax([24, 160]);          
  
    

  
  #%%
  import matplotlib.pyplot as plt;
  plt.figure(1); plt.clf();
  #plt.plot(qs);
  plt.hist(qs)
  
  
  
  
  
  #%%
  i = 20;
  i = 21;
  i = 30;
  i = 40;
  r = 25;
  center = np.unravel_index(ends[i], data.shape)
  print(center, data.shape)
  mask = binary; 
  path = con.tracePointToMask(data, mask, center, radius= r, points = special_xyz, 
                              plot = True, skel = skel, binary = binary, 
                              tubeness = None, removeLocalMask = True, maxSteps = None, verbose = False,
                              costPerDistance = 0.0);
       

  #%%

  nbs = ap.findNeighbours(ends, i, skel.shape, skel.strides, r);
  center = np.unravel_index(ends[i], skel.shape)

  nbs_xyz = np.array(np.unravel_index(nbs, skel.shape)).T
  dists = nbs_xyz - center;
  dists = np.sum(dists*dists, axis = 1);
  
  nb = np.argmin(dists)
    
  center = np.unravel_index(ends[i], data.shape)
  print(center, data.shape)
  mask = binary; 
  path = con.tracePointToNeighbor(data, mask, center, nbs_xyz[nb],
                                  radius= r, points = special_xyz, 
                                  plot = True, skel = skel, binary = binary, 
                                  tubeness = None, removeLocalMask = True, maxSteps = None, verbose = False,
                                  costPerDistance = 0.0);


  #%%

  import ClearMap.ImageProcessing.Filter.FilterKernel as fkr;
  dog = fkr.filterKernel('DoG', size = (13,13,13));
  dv.plot(dog)
  
  data_filter = ndi.correlate(np.asarray(data, dtype = float), dog);
  data_filter -= data_filter.min();
  data_filter = data_filter / 3.0;
  #dv.dualPlot(data, data_filter);
                       
  #%%add all paths
  reload(con)                              
                              
  r = 25;
  mask = binary; 
  data_new = data.copy();
  data_new[skel] = 120;
  
  skel_new = np.asarray(skel, dtype= int);
  skel_new = skel_new + binary;
  
  binary_new = binary.copy();
  
  for i,e in enumerate(special):
    center = np.unravel_index(e, data.shape);
    
    print(i, e, center)
    path   = con.tracePointToMask(data, mask, center, radius= r, points = special_xyz, 
                                  plot = False, skel = skel, binary = binary, 
                                  tubeness = None, removeLocalMask = True, maxSteps = 15000,
                                  costPerDistance = 1.0);
    
    q = con.addPathToMask(skel_new, path, value = 7);
    q = con.addPathToMask(data_new, path, value = 512);
    binary_new = con.addDilatedPathToMask(binary_new, path, iterations = 1);
 
  q = con.addPathToMask(skel_new, special_xyz, value = 6);
  for d in [3,4,5]:
    xyz = np.array(np.unravel_index(points[deg == d], data.shape)).T;
    q = con.addPathToMask(skel_new, xyz, value = d);  
  q = con.addPathToMask(data_new, special_xyz, value = 150);
  
  skel_new = skel_new + binary_new;
  try:
    con.viewer[0].setSource(skel_new);
    con.viewer[1].setSource(data_new);
  except:
    con.viewer = dv.dualPlot(skel_new, data_new);
  
  con.viewer[0].setMinMax([0, 8]);
  con.viewer[1].setMinMax([24, 160]);                
    
  
  #%%
  
  import ClearMap.ImageProcessing.Skeletonization.Skeletonize as skl
  
  skel_2 = skl.skeletonize3D(binary_new.copy());
  
  
  #%%
  
  np.save('/home/ckirst/Desktop/binarized_con.npy', binary_new)
  #%%
  
  # write image
  
  import ClearMap.IO.IO as io
  
  #r = np.asarray(128 * binary_new, dtype = 'uint8');
  #g = r.copy(); b = r.copy();
  #r[:] = r + 127 * skel_2[0];
  #g[:] = g - 128 * skel_2[0];
  #b[:] = b - 128 * skel_2[0];
  #img = np.stack((r,g,b), axis = 3)
  
  img = np.asarray(128 * binary_new, dtype = 'uint8');
  img[:] = img + 127 * skel_2[0];
  
  io.writeData('/home/ckirst/Desktop/3d.tif', img)
  
  
  #%%


  #%%

  

  
#  if removeLocalMask is not False:
#    if not removeLocalMask is True:
#      if not isinstance(removeLocalMask, tuple):
#        removeLocalMask = (removeLocalMask,);
#      removeLocalMask = removeLocalMask*3;
#      removeLocalMask = removeLocalMask[:3];
#      local_mask = extractNeighbourhood(mask_nbh, center_nbh, removeLocalMask);
#      sl = tuple([slice(c - r, c + r) for c,r in zip(center_nbh, removeLocalMask)]);
#    else:
#      local_mask = mask_nbh;
#      sl = (slice(None),) * 3;
#    
#    mask_nbh_label, mask_label = ndi.label(local_mask, structure = np.ones((3,3,3), dtype = bool));
#    ids = local_mask_label[tuple(center_nbh)] == local_mask_label;
#    
#    mask_nbh[ids] = False;
#    
#    if removeLocalTubness is not None:
#
#      
#      ker = sel.
#      ids = np.logical_and(ids, )
#      tubeness_nbh[sls] = ;