def get_orientation_motion(seg):

    brain_center = np.array(nd.center_of_mass( (seg == 5).view(np.ndarray) ),
                            dtype='float32')

    heart_center = np.array(nd.center_of_mass( (seg == 3).view(np.ndarray) ),
                                                  dtype='float32')

    left_lung = np.array(nd.center_of_mass( (seg == 1).view(np.ndarray) ),
                                              dtype='float')
    right_lung = np.array(nd.center_of_mass( (seg == 2).view(np.ndarray) ),
                                               dtype='float')

    u = brain_center - heart_center
    v = right_lung - left_lung
    
    u /= np.linalg.norm(u)

    v -= np.dot(v,u)*u
    v /= np.linalg.norm(v)

    w = np.cross(u,v)
    w /= np.linalg.norm(w)

    return ( u.astype("float32"),
             v.astype("float32"),
             w.astype("float32") )

def get_roi_center(roi_native_path, roi_mni_path):
    """Get ROI center of mass.
    Get back coordinate in img space and in coordinate space.
    Also actual center of mass.
    """
    # computations in native space
    if type(roi_native_path) is str:
        img = nib.load(roi_native_path)
    else:
        img = roi_native_path
    data = img.get_data()
    data = as_ndarray(data)
    my_map = data.copy()
    center_coords = ndimage.center_of_mass(np.abs(my_map))

    x_map, y_map, z_map = center_coords[:3]
    native_coords = np.asarray(coord_transform(x_map, y_map, z_map,
                                               img.get_affine())).tolist()
    voxel = [round(x) for x in center_coords]
    # computations in mni space
    if type(roi_mni_path) is str:
        img = nib.load(roi_mni_path)
    else:
        img = roi_mni_path
    data = img.get_data()
    data = as_ndarray(data)
    my_map = data.copy()
    mni_center_coords = ndimage.center_of_mass(np.abs(my_map))
    x_map, y_map, z_map = mni_center_coords[:3]
    mni_coords = np.asarray(coord_transform(x_map, y_map, z_map,
                                            img.get_affine())).tolist()
    # returns voxel and true center mass coords
    # returns also native and mni space coords
    return (voxel[:3], center_coords[:3], [round(x) for x in native_coords],
            [round(x) for x in mni_coords])

def guess_center_nested(image, halfwidth=50):
    '''Guess the position of the central object as two-step process

    First, this function calculates the center of mass of an image.
    This works well if the central object is the only bright source, however
    even a moderately bright source that is far away can shift the center of
    mass of an image by a few pixels. To improve the first guess the function
    selects a subimage with the halfwidth ``halfwidth`` in a second step
    and calculates the center of mass of that subimage.

    Parameters
    ----------
    image : 2d np.array
        input image
    halfwidth : int
        half width of the subimage selected in the second step.

    Returns
    -------
    xm, ym : float
        x and y coordinates estimated position of the central object
    '''
    xm, ym = ndimage.center_of_mass(np.ma.masked_invalid(image))
    n = 2 * halfwidth + 1
    subimage, xmymsmall = extract_array(image, (n, n), (xm, ym),
                                        return_position=True)
    x1, y1 = ndimage.center_of_mass(np.ma.masked_invalid(subimage))
    # xmymsmall is the xm, ym position in the coordinates of subimage
    # So, correct the initial (xm, ym) by delta(xmsmall, x1)
    return xm + (x1 - xmymsmall[0]), ym + (y1 - xmymsmall[1])

    def rebin_data(self, grid, use_psf=True):
        """Calculates the center of mass of the grid and then
        rebins so that the center pixel really is the center of the array
        For this we do a 2-d interpolation on the grid
        """
        a = psf_fitter.psffit(abs(grid), circle=False, rotate=1)
        xcen = a[2]
        ycen = a[2]
        xlen, ylen = grid.shape
        xval = arange(xlen)
        yval = arange(ylen)

        xint = interp1d(xval, self.xpos_abs)
        yint = interp1d(yval, self.ypos_abs)

        xintcen = self.xmax_pos-xint(xcen)
        yintcen = self.ymax_pos-yint(ycen)

        print self.xmax_pos, xintcen, self.ymax_pos, yintcen
        f_real = interp2d(self.xpos_rel, self.ypos_rel, real(grid))
        f_imag = interp2d(self.xpos_rel, self.ypos_rel, imag(grid))

        xnew = self.xpos_rel - xintcen
        ynew = self.ypos_rel - yintcen

        recen_grid = f_real(xnew, ynew) + 1j*f_imag(xnew, ynew)

        print nd.center_of_mass(abs(recen_grid))

        return recen_grid

    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 objectfeatures(img):
    """
    values=objectfeatures(img)

    This implements the object features described in
    "Object Type Recognition for Automated Analysis of Protein Subcellular Location"
    by Ting Zhao, Meel Velliste, Michael V. Boland, and Robert F. Murphy
    in IEEE Transaction on Image Processing
    """

    protimg = img.get("procprotein")
    dnaimg = img.channeldata.get("procdna", None)
    assert (
        dnaimg is None or protimg.shape == dnaimg.shape
    ), "pymorph.objectfeatures: DNA image is not of same size as Protein image."

    labeled, N = ndimage.label(protimg, ones((3, 3)))
    if not N:
        return np.zeros((0, 11))

    sofs = np.zeros((N, 11))
    indices = np.arange(1, N + 1)
    if dnaimg is not None:
        dnacofy, dnacofx = ndimage.center_of_mass(dnaimg)
        bindna = dnaimg > 0
        # According to the documentation, it shouldn't matter if indices is None,
        # but in my version of scipy.ndimage, you *have* to use indices.
        centers = ndimage.center_of_mass(protimg, labeled, indices)
        if N == 1:
            centers = list(centers)
        centers = np.asarray(centers)
        centers -= np.array((dnacofy, dnacofx))
        centers **= 2
        sofs[:, 1] = np.sqrt(centers.sum(1))
    locations = ndimage.find_objects(labeled, N)
    sofs[:, 9] = ndimage.measurements.sum(protimg, labeled, indices)
    for obji in xrange(N):
        slice = locations[obji]
        binobj = (labeled[slice] == (obji + 1)).copy()
        protobj = protimg[slice]
        binskel = thin(binobj)
        objhull = convexhull(binobj)
        no_of_branch_points = fast_sum(find_branch_points(binskel))
        hfeats = hullfeatures(binobj, objhull)
        sofs[obji, 0] = fast_sum(binobj)
        if dnaimg is not None:
            sofs[obji, 2] = fast_sum(binobj & bindna[slice])
        sofs[obji, 3] = hfeats[2]
        sofs[obji, 4] = euler(binobj)
        sofs[obji, 5] = hfeats[1]
        sofs[obji, 6] = fast_sum(binskel)
        sofs[obji, 7] = hfeats[0]
        sofs[obji, 9] /= fast_sum(binskel * protobj)
        sofs[obji, 10] = no_of_branch_points
    sofs[:, 2] /= sofs[:, 0]
    sofs[:, 8] = sofs[:, 6] / sofs[:, 0]
    sofs[:, 10] /= sofs[:, 6]
    return sofs

def angles2transfo(image1, image2, angleX=0, angleY=0, angleZ=0) :
    """
    Compute transformation matrix between 2 images from the angles in each directions.

    :Parameters:
     - `image1` (|SpatialImage|) -
     - `image2` (|SpatialImage|) -
     - `angleX` (int) - Rotation through angleX (degree)
     - `angleY` (int) - Rotation through angleY (degree)
     - `angleZ` (int) - Rotation through angleZ (degree)

    :Returns:
     - matrix (numpy array) - Transformation matrix
    """
    x = np.array(center_of_mass(image1))
    y = np.array(center_of_mass(image2))

    # Rx rotates the y-axis towards the z-axis
    thetaX = radians(angleX)
    Rx = np.zeros((3,3))
    Rx[0,0] = 1.
    Rx[1,1] = Rx[2,2] = cos(thetaX)
    Rx[1,2] = -sin(thetaX)
    Rx[2,1] = sin(thetaX)

    # Ry rotates the z-axis towards the x-axis
    thetaY = radians(angleY)
    Ry = np.zeros((3,3))
    Ry[0,0] = Ry[2,2] = cos(thetaY)
    Ry[0,2] = sin(thetaY)
    Ry[2,0] = -sin(thetaY)
    Ry[1,1] = 1.

    # Rz rotates the x-axis towards the y-axis
    thetaZ = radians(angleZ)
    Rz = np.zeros((3,3))
    Rz[0,0] = Rz[1,1] = cos(thetaZ)
    Rz[1,0] = sin(thetaZ)
    Rz[0,1] = -sin(thetaZ)
    Rz[2,2] = 1.

    # General rotations
    R = np.dot(np.dot(Rx,Ry),Rz)

    t = y - np.dot(R,x)

    matrix = np.zeros((4,4))
    matrix[0:3,0:3] = R
    matrix[0:3,3] = t
    matrix[2,2] = matrix[3,3] = 1.

    return matrix

def get_centers( seg ):
    brain_center = np.array(nd.center_of_mass( (seg == 2).view(np.ndarray) ),
                            dtype='float32')

    heart_center = np.array(nd.center_of_mass( (seg == 5).view(np.ndarray) ),
                                                  dtype='float32')

    left_lung = np.array(nd.center_of_mass( (seg == 3).view(np.ndarray) ),
                                              dtype='float')
    right_lung = np.array(nd.center_of_mass( (seg == 4).view(np.ndarray) ),
                                               dtype='float')

    return brain_center, heart_center, left_lung, right_lung

def align_heart(img,labels):

    BPD = get_BPD(30.0)
    CRL = get_CRL(30.0)

    brain_center = labels.ImageToWorld( np.array(nd.center_of_mass( (labels == 2).view(np.ndarray) ),
                                                  dtype='float32')[::-1] )

    heart_center =  labels.ImageToWorld( np.array(nd.center_of_mass( (labels == 5).view(np.ndarray) ),
                                                  dtype='float32')[::-1] )

    lungs_center =  labels.ImageToWorld( np.array(nd.center_of_mass(np.logical_or(labels == 3,
                                                                                   labels == 4 ).view(np.ndarray)
                                                                    ),
                                                  dtype='float32')[::-1] )

    left_lung = labels.ImageToWorld( np.array(nd.center_of_mass( (labels == 3).view(np.ndarray) ),
                                              dtype='float')[::-1] )
    right_lung = labels.ImageToWorld( np.array(nd.center_of_mass( (labels == 4).view(np.ndarray) ),
                                               dtype='float')[::-1] )
    
    u = brain_center - heart_center
    #v = lungs_center - heart_center
    v = right_lung - left_lung
    
    u /= np.linalg.norm(u)

    v -= np.dot(v,u)*u
    v /= np.linalg.norm(v)

    w = np.cross(u,v)
    w /= np.linalg.norm(w)
    
    # v = np.cross(w,u)
    # v /= np.linalg.norm(v)

    header = img.get_header()
    header['orientation'][0] = u
    header['orientation'][1] = v
    header['orientation'][2] = w

    header['origin'][:3] = heart_center
    
    header['dim'][0] = CRL
    header['dim'][1] = CRL
    header['dim'][2] = CRL
    
    new_img = img.transform( target=header, interpolation="bspline" )
    new_labels = labels.transform( target=header, interpolation="nearest" )

    return new_img, new_labels

    def com_dist(self):
        """
        This function calculates the euclidean distance between the centres
        of mass of the reference and segmentation.

        :return:
        """
        if self.flag_empty:
            return -1
        com_ref = ndimage.center_of_mass(self.ref)
        com_seg = ndimage.center_of_mass(self.seg)
        com_dist = np.sqrt(np.dot(np.square(np.asarray(com_ref) -
                                            np.asarray(com_seg)), np.square(
            self.pixdim)))
        return com_dist

def find_local_maxima(image, min_distance):
    """Find maxima in an image.

    Finds the highest-valued points in an image, such that each point is
    separted by at least min_distance.

    If there are flat regions that are all at a maxima, the enter of mass of
    the region is reported. Large flat regions of more than min_distance in
    radius will be erroneously returned as maxima even if they are not. Further
    filtering should be performed to exclude these if needed.

    Returns the position of the maxima and the value at each maximum.

    Parameters:
        image: image of arbitrary dimensionality
        min_distance: maxima found will be at least this many pixels apart

    Returns:
        centroids: list of centers of each maxima
        values: image value at each maxima

  """
    image_max = ndimage.maximum_filter(image, size=2*min_distance+1, mode='constant')
    peak_mask = (image == image_max)
    # NB: some maxima might be marked by multiple contiguous pixels if the image
    # has "plateaus". So we need to label the mask and get the centroids
    # of each of the labeled regions.
    labeled_image, num_regions = ndimage.label(peak_mask)
    label_indices = numpy.arange(1, num_regions+1)
    centroids = ndimage.center_of_mass(peak_mask, labeled_image, label_indices)
    values = ndimage.mean(image, labeled_image, label_indices)
    return numpy.array(centroids), values

def find_start_point(npa):
    print('finding start point')
#    print(npa.shape)
    len_y = npa.shape[1]-1
    j = int()
    prev = 0
    row = []
    for i in range(len_y,int(0.8*float(len_y)),-1): 
        row = npa[i,0:]
        if i<len_y:
            prev = npa[i+1,0:]
        if len(row[row>130]) and len(prev[prev>130]):
            j = i
            break
    try:
        st_pt = ndimage.center_of_mass(row)[0]
    except RuntimeWarning:
        print(row[row>130])
#    print(j,st_pt)
    try:
        pts = ndimage.measurements.center_of_mass(npa[0:int(npa.shape[1]*0.6),0:])
    except RuntimeWarning:
        print(npa[0:int(npa.shape[1]*0.6),0:])
#    print(pts)
### Testing ###
#    img = im.fromarray(npa)
#    img.convert('RGB')
#    draw = imd.Draw(img)
#    draw.ellipse((pts[1]-20,pts[0]-20,pts[1]+20,pts[0]+20),fill="red")
#    
#    img.show()
#    del draw
################
    cm_x = int(pts[1])
    return (st_pt,j,cm_x)    

def test_cmp_ndimage():
    R = (255 * np.random.rand(128, 256)).astype(np.uint16)
    R += np.arange(256)
    m0, m1 = mahotas.center_of_mass(R)
    n0, n1 = ndimage.center_of_mass(R)
    assert np.abs(n0 - m0) < 1.0
    assert np.abs(n1 - m1) < 1.0

 def calculate_life(self, cells, area, radius):
     y, x = nd.center_of_mass(cells)
     life = np.zeros(np.shape(area))
     for cell in np.transpose(area.nonzero()):
         d = np.sqrt(np.power(y - cell[0], 2) + np.power(x - cell[1], 2))
         life[cell[0], cell[1]] = (1.0 / d) * radius + random.random() * .5 + .1
     return life

def get_spec(fname, roi_start, roi_width=180, nchannels=2, force_start=False,
             **kwargs):
    """return a sipm spectrum using the cm method as a cut.
    roi_start is the star of the region of interest, + roi_width channels
    ref_cm is the reference center of mass. if None then mean(cm) of all events
    will be calculated.
    dev_cm is the allowed deviation from ref_cm. if None then std(cm)
    of all events will be calculated.
    nchannels is the number of DRS channels with data. either 1 or 2"""

    st, wd = roi_start, roi_width
    my_dtype = return_dtype(nchannels)
    if not force_start:
        st, ref_cm, dev_cm = find_start(fname, roi_start, roi_width, nchannels,
                                        **kwargs)
    else:
        cmsarr = cms_(fname, roi_start, roi_width, nchannels)
        cmhist = histogram(cmsarr, bins=512)
        ref_cm = cmhist[1][argmax(cmhist[0])]
        dev_cm = dev_cm_(cmsarr)

    with open(fname, 'r') as f:
        gen = (fromstring(event, my_dtype)[0][5]
               for event in event_generator(f, nchannels))
        specdata = [sum(event[st:st + wd]) for event in gen
                    if abs(center_of_mass(- event[st:st + wd])[0] - ref_cm)
                    < dev_cm]
    return histogram(specdata, bins=2048)

def _hull_computations(imageproc,imagehull = None):
    # Just share code between the two functions below
    if imagehull is None:
        imagehull = convexhull(imageproc > 0)

    Ahull = _bwarea(imagehull)
    Phull = _bwarea(bwperim(imagehull))

    cofy,cofx = center_of_mass(imagehull)
    hull_mu00 = imgcentmoments(imagehull,0,0,cofy,cofx)
    hull_mu11 = imgcentmoments(imagehull,1,1,cofy,cofx)
    hull_mu02 = imgcentmoments(imagehull,0,2,cofy,cofx)
    hull_mu20 = imgcentmoments(imagehull,2,0,cofy,cofx)

# Parameters of the 'image ellipse'
#   (the constant intensity ellipse with the same mass and
#   second order moments as the original image.)
#   From Prokop, RJ, and Reeves, AP.  1992. CVGIP: Graphical
#   Models and Image Processing 54(5):438-460
    hull_semimajor = sqrt((2 * (hull_mu20 + hull_mu02 + \
                    sqrt((hull_mu20 - hull_mu02)**2 + \
                    4 * hull_mu11**2)))/hull_mu00) 

    hull_semiminor = sqrt((2 * (hull_mu20 + hull_mu02 - \
                    sqrt((hull_mu20 - hull_mu02)**2 + \
                    4 * hull_mu11**2)))/hull_mu00) 
    return imagehull,Ahull, Phull, hull_semimajor, hull_semiminor

def shiftToCenter(infile,shiftfile,isEMAN=False):
	'''
	EMAN defines the rotation origin differently from other packages.
	Therefore, it needs to be recenterred according to the package
	after using EMAN proc3d rotation functions.
	'''
	# center of rotation for eman is not at length/2.
	if isEMAN:
		formatoffset = getEmanCenter()
		prefix = ''
	else:
		formatoffset = (0,0,0)
		prefix = 'non-'

	apDisplay.printMsg('Shifting map center for %sEMAN usage' % (prefix,))
	# Find center of mass of the density map
	a = mrc.read(infile)
	t = a.mean()+2*a.std()
	numpy.putmask(a,a>=t,t)
	numpy.putmask(a,a<t,0)
	center = ndimage.center_of_mass(a)
	offset = (center[0]+formatoffset[0]-a.shape[0]/2,center[1]+formatoffset[1]-a.shape[1]/2,center[2]+formatoffset[2]-a.shape[2]/2)
	offset = (-offset[0],-offset[1],-offset[2])
	apDisplay.printMsg('Shifting map center by (x,y,z)=(%.2f,%.2f,%.2f)' % (offset[2],offset[1],offset[0]))
	# shift the map
	a = mrc.read(infile)
	a = ndimage.interpolation.shift(a,offset)
	mrc.write(a,shiftfile)
	h = mrc.readHeaderFromFile(infile)
	mrc.update_file_header(shiftfile,h)

def getRealLabeledAreaCenter(image,labeled_image,indices,info):
        print "Getting real area and center"
        shape=numpy.shape(image)
        ones=numpy.ones(shape)
        area=nd.sum(ones,labels=labeled_image,index=indices)
        center=nd.center_of_mass(ones,labels=labeled_image,index=indices)

        ll=0
        try:
                len(area)
        except:
                area=[area]
                center=[center]
                try:
                        len(indices)
                except:
                        indices=[indices]
        try:
                info.keys()
        except:
                offset=1
        else:
                offset=0
        for l in indices:
                info[l-offset][0]=area[ll]
                info[l-offset][4]=center[ll]
                ll += 1
        return info     

def nucleicof(dnaimg,options=None):
    '''
    Returns a set of nuclear centres.
    '''
    labeled,N=labelnuclei(dnaimg)
    cofs=center_of_mass(dnaimg,labeled,range(1,N+1))
    return cofs

    def get_labels(self):
        """
        find clusters and extract center and size
        Parameters
        ----------

        Returns
        -------
        self.labels : 'list'
            list of cluster labels
        self.centers : 'list'
            list of cluster centers
        self.sizes : 'list'
            list of cluster sizes

        """
        b_img = self.img_b
        label_im, nb_labels = ndimage.label(b_img)

        center = np.asarray(ndimage.center_of_mass(b_img, label_im,
                            range(1, nb_labels + 1)))
        size = np.asarray(ndimage.sum(b_img, label_im,
                          range(1, nb_labels + 1)))

        self.labels = label_im
        self.centers = center
        self.sizes = size