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()

    def __iter__(self):
        """
        Return generator of (source, target, void) tuples.

        Source and target are views into a larger array. Void is a newly
        created array containing the footprint of the void.
        """
        if progress:  # pragma: no cover
            gdal.TermProgress_nocb(0)

        # analyze
        mask = (self.source == self.no_data_value)
        labels, total = ndimage.label(mask)
        items = ndimage.find_objects(labels)

        # iterate the objects
        for label, item in enumerate(items, 1):
            index = self._grow(item)       # to include the edge
            source = self.source[index]    # view into source array
            target = self.target[index]    # view into target array
            void = labels[index] == label  # the footprint of this void
            yield source, target, void

            if progress:  # pragma: no cover
                gdal.TermProgress_nocb(label / total)

def reg_median_cont(data,mask,min_size):
    emask = ndimage.morphology.binary_erosion(mask,iterations=2)
    cmask = mask[:]
    cmask[emask==1]=0
    label_im, nb_labels = ndimage.label(cmask)
    sizes = ndimage.sum(cmask,label_im,range(nb_labels+1))
    mask_size = sizes < min_size
    remove_pixel = mask_size[label_im]
    label_im[remove_pixel] = 0
    labels = np.unique(label_im)
    label_im = np.searchsorted(labels, label_im)
    labels = np.unique(label_im)
    out = np.array(label_im,dtype=np.float)
    for lab in labels:
        if( lab==0 ): continue
        try:
            slice_x, slice_y = ndimage.find_objects(label_im==lab)[0]
        except IndexError:
            print ("Bad index: "%lab)
            continue
#        print lab
        rois = data[slice_x, slice_y]
        tmask = label_im==lab 
        roim = tmask[slice_x, slice_y]
        roio = out[slice_x, slice_y]
        mean = np.ma.median(np.ma.array(rois,mask=~roim))
        roio[roim] = mean

    return out

    def make_profile_map(self, order_map, slitpos_map, lsf,
                         slitoffset_map=None):
        """
        lsf : callable object which takes (o, x, slit_pos)

        o : order (integer)
        x : detector position in dispersion direction
        slit_pos : 0..1

        x and slit_pos can be array.
        """

        iy, ix = np.indices(slitpos_map.shape)

        if slitoffset_map is not None:
            ix = ix - slitoffset_map

        profile_map = np.empty(slitpos_map.shape, "d")
        profile_map.fill(np.nan)

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

            profile1 = np.zeros(profile_map[sl].shape, "d")
            profile1[msk] = lsf(o, ix[sl][msk], slitpos_map[sl][msk])
            # TODO :make sure that renormalization is good thing to do.
            profile_sum = np.abs(profile1).sum(axis=0)
            profile1 /= profile_sum

            profile_map[sl][msk] = profile1[msk]

        return profile_map

def get_symbols(image):
  dil_eros = bin_search(dilatation_cross_numb, [image], (1, 16), 1.0, "dec")
  block_size = 50
  binary_adaptive_image = erosion(dilation(threshold_adaptive(
    array(image.convert("L")), block_size, offset=10),
      square(dil_eros)), square(dil_eros))

  all_labels = label(binary_adaptive_image, background = True)
  objects = find_objects(all_labels)

  av_width = av_height = 0
  symbols = []

  for obj in objects:
    symb = (binary_adaptive_image[obj], (obj[0].start, obj[1].start))
    symbols.append(symb)
    av_height += symb[0].shape[0]
    av_width += symb[0].shape[1]

  av_width /= float(len(objects))
  av_height /= float(len(objects))

  symbols = [symb for symb in symbols
    if symb[0].shape[0] >= av_height and symb[0].shape[1] >= av_width]

  return symbols

def DrawRectangle(img, dst, nb_labels, color=(255,0,0)):
    """
        Function definition
        +++++++++++++++++++
            
        .. py:function:: DrawRectangle(img, dst, nb_labels, color=(255,0,0))
            
            This method finds objects of interest and draws a rectangle around each one of them. 

            :param numpy_array img: image on which the objects of interest are detected.
            :param numpy_array img: image on which the rectangles are drawn.
            :param int nb_labels: number of objects to detect. Usually, the number of labels in a
                                  labeled image.
            :param tuple coler: the color of rectangles outline.
            :return: regions of interest - each one of them contains an object of interest.
            :rtype: list of lists. **Example:** rois[0][.....] contains slice_x and rois[1][.....]
                                   contains slice_y.
        """
    
    rois = []
    
    for i in range(nb_labels):
        slice_x, slice_y = ndimage.find_objects(img==i+1)[0]
        cv2.rectangle(dst,(slice_y.start,slice_x.start),(slice_y.stop,slice_x.stop),color,1)
        rois.append([slice_x, slice_y])
        
    return rois