python类mean()的实例源码

def apply(self):
        sc = self.scene
        org = self.original
        factor = self.factor

        sx, sy = sc.displacement.shape
        gx, gy = num.ogrid[0:sx, 0:sy]
        regions = sy/factor * (gx/factor) + gy/factor
        indices = num.arange(regions.max() + 1)

        def block_downsample(arr):
            res = ndimage.mean(
                arr,
                labels=regions,
                index=indices)
            res.shape = (sx/factor, sy/factor)
            return res

        sc.displacement = block_downsample(sc.displacement)
        sc.theta = block_downsample(sc.theta)
        sc.phi = block_downsample(sc.phi)
        sc.frame.dLat = org['frame.dLat'] * self.factor
        sc.frame.dLon = org['frame.dLat'] * self.factor

def downsample(voxels, step, method='max'):
    """
    downsample a voxels matrix by a factor of step. 
    downsample method options: max/mean 
    same as a pooling
    """
    assert step > 0
    assert voxels.ndim == 3 or voxels.ndim == 4
    assert method in ('max', 'mean')
    if step == 1:
        return voxels

    if voxels.ndim == 3:
        sx, sy, sz = voxels.shape[-3:]
        X, Y, Z = np.ogrid[0:sx, 0:sy, 0:sz]
        regions = sz/step * sy/step * (X/step) + sz/step * (Y/step) + Z/step
        if method == 'max':
            res = ndimage.maximum(voxels, labels=regions, index=np.arange(regions.max() + 1))
        elif method == 'mean':
            res = ndimage.mean(voxels, labels=regions, index=np.arange(regions.max() + 1))
        res.shape = (sx/step, sy/step, sz/step)
        return res
    else:
        res0 = downsample(voxels[0], step, method)
        res = np.zeros((voxels.shape[0],) + res0.shape)
        res[0] = res0
        for ind in xrange(1, voxels.shape[0]):
            res[ind] = downsample(voxels[ind], step, method)
        return res

def downsample(voxels, step, method='max'):
    """
    downsample a voxels matrix by a factor of step. 
    downsample method options: max/mean 
    same as a pooling
    """
    assert step > 0
    assert voxels.ndim == 3 or voxels.ndim == 4
    assert method in ('max', 'mean')
    if step == 1:
        return voxels

    if voxels.ndim == 3:
        sx, sy, sz = voxels.shape[-3:]
        X, Y, Z = np.ogrid[0:sx, 0:sy, 0:sz]
        regions = sz/step * sy/step * (X/step) + sz/step * (Y/step) + Z/step
        if method == 'max':
            res = ndimage.maximum(voxels, labels=regions, index=np.arange(regions.max() + 1))
        elif method == 'mean':
            res = ndimage.mean(voxels, labels=regions, index=np.arange(regions.max() + 1))
        res.shape = (sx/step, sy/step, sz/step)
        return res
    else:
        res0 = downsample(voxels[0], step, method)
        res = np.zeros((voxels.shape[0],) + res0.shape)
        res[0] = res0
        for ind in xrange(1, voxels.shape[0]):
            res[ind] = downsample(voxels[ind], step, method)
        return res

def downsample(voxels, step, method='max'):
    """
    downsample a voxels matrix by a factor of step. 
    downsample method options: max/mean 
    same as a pooling
    """
    assert step > 0
    assert voxels.ndim == 3 or voxels.ndim == 4
    assert method in ('max', 'mean')
    if step == 1:
        return voxels

    if voxels.ndim == 3:
        sx, sy, sz = voxels.shape[-3:]
        X, Y, Z = np.ogrid[0:sx, 0:sy, 0:sz]
        regions = sz/step * sy/step * (X/step) + sz/step * (Y/step) + Z/step
        if method == 'max':
            res = ndimage.maximum(voxels, labels=regions, index=np.arange(regions.max() + 1))
        elif method == 'mean':
            res = ndimage.mean(voxels, labels=regions, index=np.arange(regions.max() + 1))
        res.shape = (sx/step, sy/step, sz/step)
        return res
    else:
        res0 = downsample(voxels[0], step, method)
        res = np.zeros((voxels.shape[0],) + res0.shape)
        res[0] = res0
        for ind in xrange(1, voxels.shape[0]):
            res[ind] = downsample(voxels[ind], step, method)
        return res

def property_topomesh_area_smoothing_force(topomesh,target_areas=None):

    compute_topomesh_property(topomesh,'vertices',degree=2)
    compute_topomesh_property(topomesh,'length',degree=1)

    triangle_vertices = topomesh.wisp_property('vertices',degree=2).values()
    rotated_triangle_vertices = np.transpose([triangle_vertices[:,2],triangle_vertices[:,0],triangle_vertices[:,1]])
    antirotated_triangle_vertices = np.transpose([triangle_vertices[:,1],triangle_vertices[:,2],triangle_vertices[:,0]])
    triangle_vertices = np.append(np.append(triangle_vertices,rotated_triangle_vertices,axis=0),antirotated_triangle_vertices,axis=0)

    edge_index_list = np.array([[1, 2],[0, 1],[0, 2]])
    triangle_edge_vertices = triangle_vertices[:,edge_index_list]

    triangle_edge_vectors = topomesh.wisp_property('barycenter',degree=0).values(triangle_edge_vertices[...,1]) - topomesh.wisp_property('barycenter',degree=0).values(triangle_edge_vertices[...,0])
    #triangle_edge_lengths = np.power(np.sum(np.power(triangle_edge_vectors,2.0),axis=2),0.5)
    triangle_edge_lengths = np.linalg.norm(triangle_edge_vectors,axis=2)
    triangle_edge_directions = triangle_edge_vectors/triangle_edge_lengths[...,np.newaxis]

    triangle_perimeters = np.sum(triangle_edge_lengths,axis=1)
    triangle_areas = np.sqrt((triangle_perimeters/2.0)*(triangle_perimeters/2.0-triangle_edge_lengths[:,0])*(triangle_perimeters/2.0-triangle_edge_lengths[:,1])*(triangle_perimeters/2.0-triangle_edge_lengths[:,2]))
    triangle_areas[np.where(triangle_areas==0.0)] = 0.001

    if target_areas is None:
        target_areas = triangle_areas.mean()*np.ones(len(triangle_areas))
    else:
        target_areas = np.tile(target_areas,(3))

    area_edge_1_force = (target_areas-triangle_areas)*(triangle_edge_lengths[:,0]**2+triangle_edge_lengths[:,2]**2-triangle_edge_lengths[:,1]**2)*triangle_edge_lengths[:,1]
    area_edge_2_force = (target_areas-triangle_areas)*(triangle_edge_lengths[:,0]**2+triangle_edge_lengths[:,1]**2-triangle_edge_lengths[:,2]**2)*triangle_edge_lengths[:,2]

    area_unitary_force = -(area_edge_1_force[:,np.newaxis]*triangle_edge_directions[:,1] + area_edge_2_force[:,np.newaxis]*triangle_edge_directions[:,2])

    triangle_force = np.transpose([nd.sum(area_unitary_force[:,d],triangle_vertices[:,0],index=list(topomesh.wisps(0))) for d in xrange(3)])

    return triangle_force

def run(self, ips, snap, img, para = None):
        intenimg = WindowsManager.get(para['inten']).ips.img
        strc = ndimage.generate_binary_structure(2, 1 if para['con']=='4-connect' else 2)
        buf, n = ndimage.label(snap, strc, output=np.uint16)
        index = range(1, n+1)
        idx = (np.ones(n+1)*para['front']).astype(np.uint8)
        msk = np.ones(n, dtype=np.bool)

        if para['mean']>0: msk *= ndimage.mean(intenimg, buf, index)>=para['mean']
        if para['mean']<0: msk *= ndimage.mean(intenimg, buf, index)<-para['mean']
        if para['max']>0: msk *= ndimage.maximum(intenimg, buf, index)>=para['max']
        if para['max']<0: msk *= ndimage.maximum(intenimg, buf, index)<-para['max']
        if para['min']>0: msk *= ndimage.minimum(intenimg, buf, index)>=para['min']
        if para['min']<0: msk *= ndimage.minimum(intenimg, buf, index)<-para['min']
        if para['sum']>0: msk *= ndimage.sum(intenimg, buf, index)>=para['sum']
        if para['sum']<0: msk *= ndimage.sum(intenimg, buf, index)<-para['sum']
        if para['std']>0: msk *= ndimage.standard_deviation(intenimg, buf, index)>=para['std']
        if para['std']<0: msk *= ndimage.standard_deviation(intenimg, buf, index)<-para['std']


        xy = ndimage.center_of_mass(intenimg, buf, index)
        xy = np.array(xy).round(2).T

        idx[1:][~msk] = para['back']
        idx[0] = 0
        img[:] = idx[buf]

        WindowsManager.get(para['inten']).ips.mark = RGMark((xy.T, msk))
        WindowsManager.get(para['inten']).ips.update = True

def _roundness(x, y, pixelSize):
    """
    Implements equation (5) of Cappellari & Copin (2003)

    """
    n = x.size
    equivalentRadius = np.sqrt(n/np.pi)*pixelSize
    xBar, yBar = np.mean(x), np.mean(y)  # Geometric centroid here!
    maxDistance = np.sqrt(np.max((x - xBar)**2 + (y - yBar)**2))
    roundness = maxDistance/equivalentRadius - 1.

    return roundness

#----------------------------------------------------------------------

def _reassign_bad_bins(classe, x, y):
    """
    Implements steps (vi)-(vii) in section 5.1 of Cappellari & Copin (2003)

    """
    # Find the centroid of all successful bins.
    # CLASS = 0 are unbinned pixels which are excluded.
    #
    good = np.unique(classe[classe > 0])
    xnode = ndimage.mean(x, labels=classe, index=good)
    ynode = ndimage.mean(y, labels=classe, index=good)

    # Reassign pixels of bins with S/N < targetSN
    # to the closest centroid of a good bin
    #
    bad = classe == 0
    index = voronoi_tessellation(x[bad], y[bad], xnode, ynode, [1])
    classe[bad] = good[index]

    # Recompute all centroids of the reassigned bins.
    # These will be used as starting points for the CVT.
    #
    good = np.unique(classe)
    xnode = ndimage.mean(x, labels=classe, index=good)
    ynode = ndimage.mean(y, labels=classe, index=good)

    return xnode, ynode

#----------------------------------------------------------------------------

def property_topomesh_laplacian_epidermis_convexity_force(topomesh):
    """todo"""

    if not topomesh.has_wisp_property('vertices',degree=1,is_computed=True):
        compute_topomesh_property(topomesh,'vertices',degree=1)

    if not topomesh.has_wisp_property('barycenter',degree=3,is_computed=True):
        compute_topomesh_property(topomesh,'barycenter',degree=3)

    if not topomesh.has_wisp_property('epidermis',degree=1,is_computed=True):
        compute_topomesh_property(topomesh,'epidermis',degree=1)

    if not topomesh.has_wisp_property('cells',degree=0,is_computed=True):
        compute_topomesh_property(topomesh,'cells',degree=0)
    if not topomesh.has_wisp_property('cells',degree=1,is_computed=True):
        compute_topomesh_property(topomesh,'cells',degree=1)


    epidermis_convexity_force = array_dict(np.zeros_like(topomesh.wisp_property('barycenter',degree=0).values(),np.float32),np.array(list(topomesh.wisps(0))))

    edge_cells_degree = np.array(map(len,topomesh.wisp_property('cells',degree=1).values()))
    # edge_vertices = topomesh.wisp_property('vertices',degree=1).values()[np.where((topomesh.wisp_property('epidermis',degree=1).values()))]
    # edge_vertices = topomesh.wisp_property('vertices',degree=1).values()[np.where((topomesh.wisp_property('epidermis',degree=1).values())&(edge_cells_degree>1))]
    edge_vertices = topomesh.wisp_property('vertices',degree=1).values()[np.where((edge_cells_degree>2)|((topomesh.wisp_property('epidermis',degree=1).values())&(edge_cells_degree>1)))]

    epidermis_vertices = np.unique(edge_vertices)
    vertices_degrees = np.array(nd.sum(np.ones_like(edge_vertices),edge_vertices,index=epidermis_vertices),int)

    reversed_edge_vertices = np.transpose([edge_vertices[:,1],edge_vertices[:,0]])
    edge_vertices = np.append(edge_vertices,reversed_edge_vertices,axis=0)

    edge_vectors = topomesh.wisp_property('barycenter',degree=0).values(edge_vertices[:,1]) - topomesh.wisp_property('barycenter',degree=0).values(edge_vertices[:,0])

    laplacian_force = np.transpose([nd.sum(edge_vectors[:,0],edge_vertices[:,0],index=epidermis_vertices),
                                    nd.sum(edge_vectors[:,1],edge_vertices[:,0],index=epidermis_vertices),
                                    nd.sum(edge_vectors[:,2],edge_vertices[:,0],index=epidermis_vertices)])

    laplacian_force = laplacian_force/vertices_degrees[:,np.newaxis]

    # vertex_cell_barycenter = np.array([np.mean(topomesh.wisp_property('barycenter',degree=3).values(c),axis=0)  for c in topomesh.wisp_property('cells',0).values(epidermis_vertices)])
    # vertices_directions = topomesh.wisp_property('barycenter',degree=0).values(epidermis_vertices) - vertex_cell_barycenter
    # vertices_directions = vertices_directions/np.linalg.norm(vertices_directions,axis=1)[:,np.newaxis]
    # vertices_products = np.einsum('ij,ij->i',laplacian_force,vertices_directions)
    # convex_points = np.where(vertices_products<0.)[0]
    # laplacian_force[convex_points] -= vertices_directions[convex_points]*vertices_products[convex_points,np.newaxis]
    # laplacian_force[convex_points] = np.array([0.,0.,0.])

    epidermis_convexity_force.update(laplacian_force,keys=epidermis_vertices,erase_missing_keys=False)

    return epidermis_convexity_force.values()

def run(self, ips, imgs, para = None):
        inten = WindowsManager.get(para['inten']).ips
        if not para['slice']:
            imgs = [inten.img]
            msks = [ips.img]
        else: 
            msks = ips.imgs
            if len(msks)==1:
                msks *= len(imgs)
        buf = imgs[0].astype(np.uint16)
        strc = ndimage.generate_binary_structure(2, 1 if para['con']=='4-connect' else 2)
        idct = ['Max','Min','Mean','Variance','Standard','Sum']
        key = {'Max':'max','Min':'min','Mean':'mean',
               'Variance':'var','Standard':'std','Sum':'sum'}
        idct = [i for i in idct if para[key[i]]]
        titles = ['Slice', 'ID'][0 if para['slice'] else 1:] 
        if para['center']: titles.extend(['Center-X','Center-Y'])
        if para['extent']: titles.extend(['Min-Y','Min-X','Max-Y','Max-X'])
        titles.extend(idct)
        k = ips.unit[0]
        data, mark = [], []
        for i in range(len(imgs)):
            n = ndimage.label(msks[i], strc, output=buf)
            index = range(1, n+1)
            dt = []
            if para['slice']:dt.append([i]*n)
            dt.append(range(n))

            xy = ndimage.center_of_mass(imgs[i], buf, index)
            xy = np.array(xy).round(2).T
            if para['center']:dt.extend([xy[1]*k, xy[0]*k])

            boxs = [None] * n
            if para['extent']:
                boxs = ndimage.find_objects(buf)
                boxs = [(i[0].start, i[1].start, i[0].stop, i[1].stop) for i in boxs]
                for j in (0,1,2,3):
                    dt.append([i[j]*k for i in boxs])
            if para['max']:dt.append(ndimage.maximum(imgs[i], buf, index).round(2))
            if para['min']:dt.append(ndimage.minimum(imgs[i], buf, index).round(2))        
            if para['mean']:dt.append(ndimage.mean(imgs[i], buf, index).round(2))
            if para['var']:dt.append(ndimage.variance(imgs[i], buf, index).round(2)) 
            if para['std']:dt.append(ndimage.standard_deviation(imgs[i], buf, index).round(2))
            if para['sum']:dt.append(ndimage.sum(imgs[i], buf, index).round(2))      

            mark.append([(center, cov) for center,cov in zip(xy.T, boxs)]) 
            data.extend(list(zip(*dt)))

        IPy.table(inten.title+'-region statistic', data, titles)
        inten.mark = Mark(mark)
        inten.update = True

def _cvt_equal_mass(x, y, signal, noise, xnode, ynode, pixelsize, quiet, sn_func, wvt):
    """
    Implements the modified Lloyd algorithm
    in section 4.1 of Cappellari & Copin (2003).

    NB: When the keyword WVT is set this routine includes
    the modification proposed by Diehl & Statler (2006).

    """
    dens2 = (signal/noise)**4     # See beginning of section 4.1 of CC03
    scale = np.ones_like(xnode)   # Start with the same scale length for all bins

    for it in range(1, xnode.size):  # Do at most xnode.size iterations

        xnode_old, ynode_old = xnode.copy(), ynode.copy()
        classe = voronoi_tessellation(x, y, xnode, ynode, scale)

        # Computes centroids of the bins, weighted by dens**2.
        # Exponent 2 on the density produces equal-mass Voronoi bins.
        # The geometric centroids are computed if WVT keyword is set.
        #
        good = np.unique(classe)
        if wvt:
            for k in good:
                index = np.flatnonzero(classe == k)   # Find subscripts of pixels in bin k.
                xnode[k], ynode[k] = np.mean(x[index]), np.mean(y[index])
                sn = sn_func(index, signal, noise)
                scale[k] = np.sqrt(index.size/sn)  # Eq. (4) of Diehl & Statler (2006)
        else:
            mass = ndimage.sum(dens2, labels=classe, index=good)
            xnode = ndimage.sum(x*dens2, labels=classe, index=good)/mass
            ynode = ndimage.sum(y*dens2, labels=classe, index=good)/mass

        diff2 = np.sum((xnode - xnode_old)**2 + (ynode - ynode_old)**2)
        diff = np.sqrt(diff2)/pixelsize

        if not quiet:
            print('Iter: %4i  Diff: %.4g' % (it, diff))

        if diff < 0.1:
            break

    # If coordinates have changed, re-compute (Weighted) Voronoi Tessellation of the pixels grid
    #
    if diff > 0:
        classe = voronoi_tessellation(x, y, xnode, ynode, scale)
        good = np.unique(classe)  # Check for zero-size Voronoi bins

    # Only return the generators and scales of the nonzero Voronoi bins

    return xnode[good], ynode[good], scale[good], it

#-----------------------------------------------------------------------