Python ndimage.find_objects函数代码示例

def align_image_to_ellipse(coeffs, image):
    """
    Given the coefficients of an ellipse in 2D and a binary 
    image, return the angle required to align the image to the
    principal axes of the ellipse (with the longest axis
    as the first major 'hump' on the left).
    """
    
    coeff_a, coeff_b, coeff_c = coeffs[:3]
        
    # Calculate tan(angle) for the angle of rotation of the major axis
    preangle = coeff_b / (coeff_a - coeff_c)
    
    if not np.isinf(preangle):
        # Take the arctan and convert to degrees, which is what 
        # ndimage.rotate uses.
        angle = radians_to_degrees(-0.5 * np.arctan(preangle))
        
        # Order = 0 prevents interpolation from being done and screwing 
        # with our object boundaries.
        rotated = ndimage.rotate(image, angle, order=0)
        
        # Pull out the height/width of just the object.
        try:    
            height, width = rotated[ndimage.find_objects(rotated)[0]].shape
        except IndexError:
            raise EllipseAlignmentError("Can't find object after " \
                + "initial rotation.")
    else:
        angle = 0.
        height, width = image.shape
    
    # we want the height (first axis) to be the major axis.
    if width > height:
        angle -= 90.0
        rotated = ndimage.rotate(image, angle, order=0)
    
    # Correct so that in budding cells, the "major" hump is always
    # on the first.          
    if np.argmax(rotated.sum(axis=1)) > rotated.shape[0] // 2:
        angle -= 180.0
        rotated = ndimage.rotate(image, angle, order=0)
    
    # Do a find_objects on the resultant array after rotation in 
    # order to _just_ get the object and not any of the extra 
    # space that's been added.
    try:
        bounds = ndimage.find_objects(rotated)[0]
    except IndexError:
        raise EllipseAlignmentError("Can't find object after final rotation.")
    
    return rotated[bounds], angle

  def find_local_maxima(self, data, neighborhood_size):
    """ 
     find local maxima within neighborhood 
      idea from http://stackoverflow.com/questions/9111711
      (get-coordinates-of-local-maxima-in-2d-array-above-certain-value)
    """

    # find local maxima in image (width specified by neighborhood_size)
    data_max = filters.maximum_filter(data,neighborhood_size);
    maxima   = (data == data_max);
    assert np.sum(maxima) > 0;        # we should always find local maxima
  
    # remove connected pixels (plateaus)
    labeled, num_objects = ndimage.label(maxima)
    slices = ndimage.find_objects(labeled)
    maxima *= 0;
    for dx,dy in slices:
      maxima[(dx.start+dx.stop-1)/2, (dy.start+dy.stop-1)/2] = 1

    # calculate difference between local maxima and lowest 
    # pixel in neighborhood (will be used in select_local_maxima)
    data_min = filters.minimum_filter(data,neighborhood_size);
    diff     = data_max - data_min;
    self._maxima = maxima;
    self._diff   = diff;

    return maxima,diff

    def find_albino_features(self, T, im):
        import scipy.ndimage as ndi

        binarized = zeros_like(T)
        binarized[T > self.albino_threshold] = True
        (labels, nlabels) = ndi.label(binarized)
        slices = ndi.find_objects(labels)

        intensities = []
        transform_means = []

        if len(slices) < 2:
            return (None, None)

        for s in slices:

            transform_means.append(mean(T[s]))
            intensities.append(mean(im[s]))

        sorted_transform_means = argsort(transform_means)
        candidate1 = sorted_transform_means[-1]
        candidate2 = sorted_transform_means[-2]

        c1_center = array(ndi.center_of_mass(im, labels, candidate1 + 1))
        c2_center = array(ndi.center_of_mass(im, labels, candidate2 + 1))

        if intensities[candidate1] > intensities[candidate2]:
            return (c2_center, c1_center)
        else:
            return (c1_center, c2_center)

def precomputestats(image):
    image.lazy_load()
    image.temp['bgsubprotein'] = bgsub(image.channeldata['protein'].copy())
    if 'dna' in image.channeldata:
        image.temp['bgsubdna'] = bgsub(image.channeldata['dna'].copy())
    if image.regions is not None:
        image.temp['region_ids'] = ndimage.find_objects(image.regions)

    def edges(cls):
        from scipy import ndimage, misc
        import numpy as np
        from skimage import feature
        col = Image.open("f990.jpg")
        gray = col.convert('L')

        # Let numpy do the heavy lifting for converting pixels to pure black or white
        bw = np.asarray(gray).copy()

        # Pixel range is 0...255, 256/2 = 128
        bw[bw < 245]  = 0    # Black
        bw[bw >= 245] = 255 # White
        bw[bw == 0] = 254
        bw[bw == 255] = 0
        im = bw
        im = ndimage.gaussian_filter(im, 1)
        edges2 = feature.canny(im, sigma=2)
        labels, numobjects =ndimage.label(im)
        slices = ndimage.find_objects(labels)
        print('\n'.join(map(str, slices)))
        misc.imsave('f990_sob.jpg', im)
        return

        #im = misc.imread('f990.jpg')
        #im = ndimage.gaussian_filter(im, 8)
        sx = ndimage.sobel(im, axis=0, mode='constant')
        sy = ndimage.sobel(im, axis=1, mode='constant')
        sob = np.hypot(sx, sy)
        misc.imsave('f990_sob.jpg', edges2)

    def extract_slit_profile(self, order_map, slitpos_map, data,
                             x1, x2, bins=None):

        x1, x2 = int(x1), int(x2)

        slices = ni.find_objects(order_map)
        slit_profile_list = []
        if bins is None:
            bins = np.linspace(0., 1., 40)

        for o in self.orders:
            sl = slices[o-1][0], slice(x1, x2)
            msk = (order_map[sl] == o)

            #ss = slitpos_map[sl].copy()
            #ss[~msk] = np.nan

            d = data[sl][msk]
            finite_mask = np.isfinite(d)
            hh = np.histogram(slitpos_map[sl][msk][finite_mask],
                              weights=d[finite_mask], bins=bins,
                              )
            slit_profile_list.append(hh[0])

        return bins, slit_profile_list

 def segment_with_label(self, img):
     # Next-nearest neighbors
     struct_nnn = np.ones((3, 3), dtype=int)
     labels, _ = ndimage.label(img, structure=struct_nnn)
     # np.savetxt(c.temp_path('labels.txt'), labels, fmt='%d')
     object_slices = ndimage.find_objects(labels)
     return labels, object_slices

def label_components(img):
    # label connected components in image.
    labeled_img, count = spimg.label(img)
    # obtain list of tuple, which are slices of array with distinct label
    slices = spimg.find_objects(labeled_img)

    return labeled_img, slices

def myfindChessboardCorners(im,dim):
    gr=30
    patern=np.zeros((gr,gr),dtype='uint8')
    patern[:gr/2,:gr/2]=255
    patern[gr/2:,gr/2:]=255
    m1=cv2.matchTemplate(im,patern,cv2.TM_CCORR_NORMED)
    patern=np.ones((gr,gr),dtype='uint8')*255
    patern[:gr/2,:gr/2]=0
    patern[gr/2:,gr/2:]=0
    m2=cv2.matchTemplate(im,patern,cv2.TM_CCORR_NORMED)
    #m=np.bitwise_or(m1>0.9,m2>0.9)
    #import pdb;pdb.set_trace()
    tresh=0.95
    labels=ndimage.label(np.bitwise_or(m1>tresh,m2>tresh))
    if labels[1]!=dim[0]*dim[1]:
        return False,[]
    objs=ndimage.find_objects(labels[0])
    corners=[]
    for xx,yy in objs:
        xpos=(xx.start+xx.stop)/2.0#+gr/2-0.5
        ypos=(yy.start+yy.stop)/2.0#+gr/2-0.5
        se=5
        #import pdb;pdb.set_trace()
        minVal, maxVal, minLoc, maxLoc=cv2.minMaxLoc(m2[xpos-se:xpos+se,ypos-se:ypos+se])
        if maxVal<tresh:
            minVal, maxVal, minLoc, maxLoc=cv2.minMaxLoc(m1[xpos-se:xpos+se,ypos-se:ypos+se])
        xpos+=-se+maxLoc[0]+gr/2-0.5
        ypos+=-se+maxLoc[1]+gr/2-0.5
        
        #xpos=xx.start+gr/2
        #ypos=yy.start+gr/2
        corners.append((ypos,xpos) )
    return True,np.array(corners)

def get_stomata(max_proj_image, min_obj_size=200, max_obj_size=1000):
    """Performs image segmentation from a max_proj_image.
     Disposes of objects in range min_obj_size to
    max_obj_size

    :param max_proj_image: the maximum projection image
    :type max_proj_image: numpy.ndarray, uint16
    :param min_obj_size: minimum size of object to keep
    :type min_obj_size: int
    :param max_obj_size: maximum size of object to keep
    :type max_obj_size: int
    :returns: list of [ [coordinates of kept objects - list of slice objects],
                        binary object image - numpy.ndarray,
                        labelled object image - numpy.ndarray
                     ]

    """

    # pore_margin = 10
    # max_obj_size = 1000
    # min_obj_size = 200
    # for prop, value in segment_options:
    #     if prop == 'pore_margin':
    #         pore_margin = value
    #     if prop == 'max_obj_size':
    #         max_obj_size = value
    #     if prop == 'min_obj_size':
    #         min_obj_size = value
    #
    # print(pore_margin)
    # print(max_obj_size)
    # print(min_obj_size)

    #rescale_min = 50
    #rescale_max= 100
    #rescaled = exposure.rescale_intensity(max_proj_image, in_range=(rescale_min,rescale_max))
    rescaled = max_proj_image
    seed = np.copy(rescaled)
    seed[1:-1, 1:-1] = rescaled.max()
    #mask = rescaled
    #if gamma != None:
    #    rescaled = exposure.adjust_gamma(max_proj_image, gamma)
    #filled = reconstruction(seed, mask, method='erosion')
    closed = dilation(rescaled)
    seed = np.copy(closed)
    seed[1:-1, 1:-1] = closed.max()
    mask = closed


    filled = reconstruction(seed, mask, method='erosion')
    label_objects, nb_labels = ndimage.label(filled)
    sizes = np.bincount(label_objects.ravel())
    mask_sizes = sizes
    mask_sizes = (sizes > min_obj_size) & (sizes < max_obj_size)
    #mask_sizes = (sizes > 200) & (sizes < 1000)
    mask_sizes[0] = 0
    big_objs = mask_sizes[label_objects]
    stomata, _ = ndimage.label(big_objs)
    obj_slices = ndimage.find_objects(stomata)
    return [obj_slices, big_objs, stomata]

 def autofocus(self, data, zsp, zsu, zind, inclusionList,
               inclusionDict, lastLabelData):
     labelData = ndimage.label(zsp.binary[zind], output=numpy.int32)[0]
     slicess = ndimage.find_objects(labelData)
     for label, slices in enumerate(slicess, start=1):
         slices = [slice(max(sli.start - 4, 0), sli.stop + 4) \
                   for sli in slices]
         dataDetail = data[slices].astype(float)
         footprint = labelData[slices] == label
         newInclusion = Inclusion(zind, label, slices, dataDetail,
                                  footprint, zsp, zsu)
         if not newInclusion.valid:
             labelData[slices] = numpy.where(
                 footprint, 0, labelData[slices]
                 )
             continue
         inclusionList.append(newInclusion)
         inclusionDict[newInclusion.index] = newInclusion
         if lastLabelData is None:
             continue
         lastLabelDetail = numpy.where(footprint,
                                       lastLabelData[slices], 0)
         for l in numpy.unique(lastLabelDetail)[1:]:
             inclusionDict[(zind - 1, l)].append(newInclusion)
     return labelData

def short_branches():
    """
    Visualization of short branches of the skeleton.
    
    """
    data1_sk = glob.glob('/backup/yuliya/vsi05/skeletons_largdom/*.h5')
    data1_sk.sort()

    for i,j, k in zip(d[1][37:47], data1_sk[46:56], ell[1][37:47]):
        g = nx.read_gpickle(i)
        dat = tb.openFile(j)
        skel = np.copy(dat.root.skel)
        bra = np.copy(dat.root.branches)
        mask = np.zeros_like(skel)    
        dat.close()
    
        length = nx.get_edge_attributes(g, 'length')
        number = nx.get_edge_attributes(g, 'number')
        num_dict = {}
        for m in number:
            for v in number[m]:
                num_dict.setdefault(v, []).append(m)
        find_br = ndimage.find_objects(bra)
        for l in list(length.keys()):
            if length[l]<0.5*k: #Criteria
                for b in number[l]:
                    mask[find_br[b-1]] = bra[find_br[b-1]]==b
        mlab.figure(bgcolor=(1,1,1), size=(1200,1200))
        mlab.contour3d(skel, colormap='hot')
        mlab.contour3d(mask)
        mlab.savefig('/backup/yuliya/vsi05/skeletons/short_bran/'+ i[42:-10] + '.png')
        mlab.close()

    def RetrieveObjects(self, masterChan=0, orientChan=1, orient_dir=-1):
        objs = ndimage.find_objects(self.image.labels)

        self.objects = [
            BlobObject(self.GetRegion(i, objs), masterChan, orientChan, orient_dir)
            for i in range(self.image.labels.max())
        ]

    def GetRegion(self, index, objects=None):
        if not objects:
            objects = ndimage.find_objects(self.image.labels)

        o = objects[index]

        mask = self.image.labels[o] == (index + 1)

        slx, sly, slz = o

        X, Y, Z = np.ogrid[slx, sly, slz]
        vs = (
            1e3 * self.image.mdh["voxelsize.x"],
            1e3 * self.image.mdh["voxelsize.y"],
            1e3 * self.image.mdh["voxelsize.z"],
        )

        return [
            DataBlock(
                np.maximum(self.image.data[slx, sly, slz, j] - self.image.data[slx, sly, slz, j].min(), 0) * mask,
                X,
                Y,
                Z,
                vs,
            )
            for j in range(self.image.data.shape[3])
        ]

    def sum_peaks(self, width):
        """
        Find peaks, then sum area around them for whole stack.

        If we're going to do this _properly_ we need a way to find areas that
        _don't_ have any beads nearby inorder to calculate noise and offset.
        """
        # fit the blobs first to find valid spots
        my_peaks = self.peakfinder

        peakfits = my_peaks.fit_blobs(diameter=width)
        # now reset the blobs to the fit values
        my_peaks.blobs = peakfits[['y0', 'x0', 'sigma_x', 'amp']].values

        # label again
        my_labels = my_peaks.label_blobs(diameter=width)

        # find all the objects.
        my_objects = ndi.find_objects(my_labels)

        my_medians = np.median(self.data, axis=(1, 2))

        my_sums = np.array([self.data[:, obj[0], obj[1]].sum((1, 2))
                            for obj in my_objects])

        self.sums = my_sums - my_medians
        # reset blobs to original
        self.peakfinder.find_blobs()