我们从Python开源项目中，提取了以下7个代码示例，用于说明如何使用scipy.ndimage.center_of_mass()。
def geometric_L1(spia): """ """ background = spia._background L1_labels = spia.cell_first_layer() L1_cells_bary = spia.center_of_mass(L1_labels, verbose=True) background_neighbors = spia.neighbors(L1_labels, min_contact_area=10., real_area=True) background_neighbors = set(background_neighbors) & set(L1_labels) L1_cells_bboxes = spia.boundingbox(L1_labels) print "-- Searching for the median voxel of each epidermis wall ..." dict_wall_voxels, epidermis_wall_median, median2bary_dist = {}, {}, {} for label_2 in background_neighbors: dict_wall_voxels[(background,label_2)] = wall_voxels_between_two_cells(spia.image, background, label_2, bbox = L1_cells_bboxes[label_2], verbose = False) epidermis_wall_median[label_2] = find_wall_median_voxel(dict_wall_voxels[(background,label_2)], verbose = False) median2bary_dist[label_2] = distance(L1_cells_bary[label_2], epidermis_wall_median[label_2]) return median2bary_dist, epidermis_wall_median, L1_cells_bary
def find_peaks_minmax(z, separation=5., threshold=10.): """ Method to locate the positive peaks in an image by comparing maximum and minimum filtered images. Parameters ---------- z : ndarray Matrix of image intensities. separation : float Expected distance between peaks. threshold : float Minimum difference between maximum and minimum filtered images. Returns ------- peaks: array with dimensions (npeaks, 2) that contains the x, y coordinates for each peak found in the image. """ data_max = ndi.filters.maximum_filter(z, separation) maxima = (z == data_max) data_min = ndi.filters.minimum_filter(z, separation) diff = ((data_max - data_min) > threshold) maxima[diff == 0] = 0 labeled, num_objects = ndi.label(maxima) peaks = np.array( ndi.center_of_mass(z, labeled, range(1, num_objects + 1))) return clean_peaks(peaks)
def autofind_pois(self, neighborhood_size=1, min_threshold=10000, max_threshold=1e6): """Automatically search the xy scan image for POIs. @param neighborhood_size: size in microns. Only the brightest POI per neighborhood will be found. @param min_threshold: POIs must have c/s above this threshold. @param max_threshold: POIs must have c/s below this threshold. """ # Calculate the neighborhood size in pixels from the image range and resolution x_range_microns = np.max(self.roi_map_data[:, :, 0]) - np.min(self.roi_map_data[:, :, 0]) y_range_microns = np.max(self.roi_map_data[:, :, 1]) - np.min(self.roi_map_data[:, :, 1]) y_pixels = len(self.roi_map_data) x_pixels = len(self.roi_map_data[1, :]) pixels_per_micron = np.max([x_pixels, y_pixels]) / np.max([x_range_microns, y_range_microns]) # The neighborhood in pixels is nbhd_size * pixels_per_um, but it must be 1 or greater neighborhood_pix = int(np.max([math.ceil(pixels_per_micron * neighborhood_size), 1])) data = self.roi_map_data[:, :, 3] data_max = filters.maximum_filter(data, neighborhood_pix) maxima = (data == data_max) data_min = filters.minimum_filter(data, 3 * neighborhood_pix) diff = ((data_max - data_min) > min_threshold) maxima[diff is False] = 0 labeled, num_objects = ndimage.label(maxima) xy = np.array(ndimage.center_of_mass(data, labeled, range(1, num_objects + 1))) for count, pix_pos in enumerate(xy): poi_pos = self.roi_map_data[pix_pos[0], pix_pos[1], :][0:3] this_poi_key = self.add_poi(position=poi_pos, emit_change=False) self.rename_poi(poikey=this_poi_key, name='spot' + str(count), emit_change=False) # Now that all the POIs are created, emit the signal for other things (ie gui) to update self.signal_poi_updated.emit()
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 inertia_axis(self, labels = None, real = True, verbose=False): """ Return the inertia axis of cells, also called the shape main axis. Return 3 (3D-oriented) vectors by rows and 3 (length) values. """ # Check the provided `labels`: labels = self.label_request(labels) # results inertia_eig_vec = [] inertia_eig_val = [] N = len(labels); percent=0 for i,label in enumerate(labels): if verbose and i*100/float(N) >= percent: print "{}%...".format(percent),; percent += 10 if verbose and i+1==N: print "100%" slices = self.boundingbox(label, real=False) center = copy.copy(self.center_of_mass(label, real=False)) # project center into the slices sub_image coordinate if slices is not None: for i,slice in enumerate(slices): center[i] = center[i] - slice.start label_image = (self.image[slices] == label) else: print 'No boundingbox found for label {}'.format(label) label_image = (self.image == label) # compute the indices of voxel with adequate label xyz = label_image.nonzero() if len(xyz)==0: continue # obviously no reasons to go further ! coord = coordinates_centering3D(xyz, center) # compute the variance-covariance matrix (1/N*P.P^T): cov = compute_covariance_matrix(coord) # Find the eigen values and vectors. eig_val, eig_vec = eigen_values_vectors(cov) # convert to real-world units if asked: if real: for i in xrange(3): eig_val[i] *= np.linalg.norm( np.multiply(eig_vec[i],self._voxelsize) ) inertia_eig_vec.append(eig_vec) inertia_eig_val.append(eig_val) if len(labels)==1 : return return_list_of_vectors(inertia_eig_vec[0]), inertia_eig_val[0] else: return self.convert_return(return_list_of_vectors(inertia_eig_vec),labels), self.convert_return(inertia_eig_val,labels)
def reduced_inertia_axis(self, labels = None, real = True, verbose=False): """ Return the REDUCED (centered coordinates standardized) inertia axis of cells, also called the shape main axis. Return 3 (3D-oriented) vectors by rows and 3 (length) values. """ # Check the provided `labels`: labels = self.label_request(labels) # results inertia_eig_vec = [] inertia_eig_val = [] N = len(labels); percent=0 for i,label in enumerate(labels): if verbose and i*100/float(N) >= percent: print "{}%...".format(percent),; percent += 10 if verbose and i+1==N: print "100%" slices = self.boundingbox(label, real=False) center = copy.copy(self.center_of_mass(label, real=False)) # project center into the slices sub_image coordinate if slices is not None: for i,slice in enumerate(slices): center[i] = center[i] - slice.start label_image = (self.image[slices] == label) else: print 'No boundingbox found for label {}'.format(label) label_image = (self.image == label) # compute the indices of voxel with adequate label xyz = label_image.nonzero() if len(xyz)==0: continue # obviously no reasons to go further ! coord = coordinates_centering3D(xyz, center) # compute the variance-covariance matrix (1/N*P.P^T): cov = compute_covariance_matrix(coord) # Find the eigen values and vectors. eig_val, eig_vec = eigen_values_vectors(cov) # convert to real-world units if asked: if real: for i in xrange(3): eig_val[i] *= np.linalg.norm( np.multiply(eig_vec[i],self._voxelsize) ) inertia_eig_vec.append(eig_vec) inertia_eig_val.append(eig_val) if len(labels)==1 : return return_list_of_vectors(inertia_eig_vec[0],by_row=1), inertia_eig_val[0] else: return self.convert_return(return_list_of_vectors(inertia_eig_vec,by_row=1),labels), self.convert_return(inertia_eig_val,labels)
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
