individual_N2_X_Y_cluster_stage_watershed.py 文件源码

def makeSegmentation(d, s, nbins = 256, verbose = True, sigma = 0.75, min_distance = 1):
  fn = os.path.join(datadir, 'data_stage%d_%s.npy' % (s,features[0]));
  dists = np.load(fn);
  fn = os.path.join(datadir, 'data_stage%d_%s.npy' % (s,features[1]));
  rots = np.load(fn);

  ddata = np.vstack([np.log(dists[:,d]), (rots[:,d])]).T
  #gmmdata = np.vstack([dists[:,j], (rots[:,j])]).T
  #ddata.shape

  nanids = np.logical_or(np.any(np.isnan(ddata), axis=1), np.any(np.isinf(ddata), axis=1));
  ddata = ddata[~nanids,:];
  #ddata.shape
  imgbin = None;
  img2 = smooth(ddata, nbins = [nbins, nbins], sigma = (sigma,sigma))
  #img = smooth(ddata, nbins = [nbins, nbins], sigma = (1,1))

  local_maxi = peak_local_max(img2, indices=False, min_distance = min_distance)
  imgm2 = img2.copy();
  imgm2[local_maxi] = 3 * imgm2.max();

  if verbose:
    imgbin = smooth(ddata, nbins = [nbins, nbins], sigma = None)
    plt.figure(220); plt.clf()
    plt.subplot(2,2,1)
    plt.imshow(imgbin)
    plt.subplot(2,2,2)
    plt.imshow(img2)
    plt.subplot(2,2,3);
    plt.imshow(imgm2,  cmap=plt.cm.jet, interpolation='nearest')

  markers = ndi.label(local_maxi)[0]
  labels = watershed(-img2, markers, mask = None);
  print "max labels: %d" % labels.max()

  if verbose:
    fig, axes = plt.subplots(ncols=3, sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'})
    ax0, ax1, ax2 = axes
    ax0.imshow(img2, cmap=plt.cm.jet, interpolation='nearest')
    ax0.set_title('PDF')
    labels[imgbin==0] = 0;
    labels[0,0] = -1;
    ax1.imshow(labels, cmap=plt.cm.rainbow, interpolation='nearest')
    ax1.set_title('Segmentation on Data')
    #labelsws[0,0] = -1;
    #ax2.imshow(labelsws, cmap=plt.cm.rainbow, interpolation='nearest')
    #ax1.set_title('Segmentation Full')

  return labels;


#classification for a specific work based on the segmentation above