Scripts: volumecenter.py

#
# Place a marker at the center of mass of a map for the region above the
# displayed contour level.
#
# Open a map file and then open this script to place the center marker.
#
# To place markers on more than one map, choose the map in the volume dialog
# by clicking on its histogram or name (so the name above the histogram is
# highlighted) and then open this script (File / Open...).
#

# Get highlighted volume from volume dialog
import VolumeViewer as VV
v = VV.active_volume()

# Find highest displayed contour level.
level = max(v.surface_levels)

# Get 3-d array of map values.
m = v.data.full_matrix()

# Find indices of map values above displayed threshold.
kji = (m >= level).nonzero()

# Compute total mass above threshold.
msum = m[kji].sum()

# Compute mass-weighted center
center = [(i*m[kji]).sum()/msum for i in kji]
center.reverse()        # k,j,i -> i,j,k index order
print 'Center of mass grid index for', v.name, center

# Convert index units to physical units.
c = v.data.ijk_to_xyz(center)

# Create marker model with same id number as volume.
import VolumePath as VP
mset = VP.Marker_Set(v.name + ' center')
mset.marker_molecule(model_id = (v.id, v.subid))

# Place marker at center position.
color = (1,1,0,1)       # (red, green, blue, opacity).
mset.place_marker(c, rgba = color, radius = min(v.data.step))