dots.py

   1
   2 import pandas as pd
   3 import numpy as np
   4 import libtiff
   5 import skimage.exposure
   6 import skimage.feature
   7 import skimage.filters
   8 import skimage.morphology
   9
  10 def make_mask(zstack, nbins=32768):
  11     thresh = skimage.filters.threshold_otsu(zstack, nbins=32768)
  12     return zstack > thresh
  13
  14 def build_blob_array(stack):
  15     data = [[[False for i in range(len(stack[0][0]))] for j in range(len(stack[0]))] for k in range(len(stack))]
  16     for im in range(len(stack))): #was stack
  17         print "processing " + str(im)
  18         yx = skimage.exposure.equalize_adapthist(stack[im], clip_limit=0.00025, nbins=32768)
  19
  20         zmx = np.max(yx)
  21         zmn = np.min(yx)
  22         yx = (yx - zmn) / (zmx - zmn)
  23
  24         blobs_log = skimage.feature.blob_log(yx, min_sigma=1, max_sigma=3, num_sigma=20, threshold=0.025)
  25         print str(len(blobs_log)) + " blobs found"
  26         if len(blobs_log) > 1:
  27             blobs_log[:, 2] = blobs_log[:, 2] * math.sqrt(2)
  28             for i in blobs_log:
  29                 data[im][int(i[0])][int(i[1])] = True
  30     return data
  31
  32 def make_watershed(mask):
  33     distance = ndi.distance_transform_edt(mask)
  34     local_max = skimage.feature.peak_local_max(distance, indices=False, footprint=np.ones((31, 31)), labels=mask)
  35     markers = ndi.label(local_max)[0]
  36     labels = skimage.morphology.watershed(-distance, markers, mask=mask)
  37     return labels
  38
  39 def make_gene_layers(gene_tif):
  40     gene_z = gene_tif.get_tiff_array()
  41     gene_s = np.arange(0, len(gene_z)+1, len(gene_z) / 6)
  42     gene_r = zip(gene_s[:-1], (gene_s + 1)[1:])
  43     gene_y = np.array(gene_z[gene_r[2][0]:gene_r[2][1]], dtype=np.uint16)
  44     return gene_y
  45
  46 def make_nissl_watershed(nissl_tif):
  47     nissl_z = nissl_tif.get_tiff_array()
  48     nissl_s = np.arange(0, len(nissl_z)+1, len(nissl_z) / 6)
  49     nissl_r = zip(nissl_s[:-1], (nissl_s + 1)[1:])
  50     nissl_y = np.array(nissl_z[nissl_r[1][0]:nissl_r[1][1]], dtype=np.uint16)
  51     nissl_watershed = [[] for k in range(len(nissl_y))]
  52
  53     for i in range(len(nissl_y)):
  54         mask = make_mask(nissl_y[i])
  55         watershed = make_watershed(mask)
  56         nissl_watershed[i] = watershed
  57     return nissl_watershed
  58
  59 def correct_labels(nissl_watershed):
  60     nissl_watershed_corrected = [[[0 for i in range(len(nissl_watershed[0][0]))] for j in range(len(nissl_watershed[0]))] for k in range(len(nissl_watershed))]
  61     nissl_watershed_corrected[0] = [i for i in nissl_watershed[0]]
  62     nissl_watershed_corrected = np.asarray(nissl_watershed_corrected)
  63
  64     for k in range(1, len(nissl_watershed_corrected)):
  65         cell_number_max = nissl_watershed_corrected[k-1].max() #changed to corrected to keep accurate count
  66         cell_number_two_max = nissl_watershed[k].max()
  67         and_map = np.logical_and(nissl_watershed[k-1], nissl_watershed[k])
  68         and_map_pos = np.where(and_map)
  69         cell_set = set(nissl_watershed[k-1][np.where(and_map)])
  70         inter_layer_dict = dict()
  71         for i in range(len(and_map_pos[0])):
  72             inter_layer_dict[nissl_watershed[k][and_map_pos[0][i]][and_map_pos[1][i]]] = nissl_watershed[k-1][and_map_pos[0][i]][and_map_pos[1][i]]
  73         counter = 0
  74         inter_layer_dict[0] = 0
  75         for i in range(1, cell_number_two_max + 1):
  76             if not i in inter_layer_dict:
  77                 counter += 1
  78                 inter_layer_dict[i] = cell_number_max + counter
  79
  80
  81         for j in range(1024):
  82             for i in range(1024):
  83                 nissl_watershed_corrected[k][i][j] = inter_layer_dict[nissl_watershed[k][i][j]]
  84
  85     return nissl_watershed_corrected