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
python类center_of_mass()的实例源码
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
