Scripts: blend.py

#
# Blend two or more solid volume renderings creating a new volumetric model.
#
#  runscript blend.py #1-3 #5
#
# Maps must be shown in solid style and be the same size, and the color red, green and blue components
# and opacity values are added, not averaged.  The first argument lists the maps to blend, and the second
# argument is the id number of the resulting model.  The result is not a map.  The result will not update
# if the map display is changed.  To update, redisplay the maps as desired then rerun this script.
#
# Script made for James Fraser to blend two x-ray maps in complementary green/magenta colors
# so that overlap regions appear gray, and differences appear green or magenta.
# An example is figure 6 in this article
#
#   http://www.sciencedirect.com/science/article/pii/S1047847704001236
#
def blend_maps(vlist, vresult):

    # Position result to map first blended map.
    s0 = vlist[0].solid
    vresult.set_array_coordinates(s0.transform)

    # Get color plane callbacks.
    planes = []
    for v in vlist:
        s = v.solid
        s.color_mode = s.c_mode = 'rgba8'
        cmap, cmap_range = s.colormap()
        planes.append(lambda axis, plane, s=s, cmap=cmap, cmap_range=cmap_range:
                      s.color_values(axis, plane, cmap, cmap_range))

    def get_color_plane(axis, plane, planes = planes):
        p = None
        for pfunc in planes:
            vp = pfunc(axis, plane)
            if p is None:
                p = vp.copy()
            else:
                p += vp
        return p
    vresult.set_color_plane_callback(s0.size, get_color_plane)

def blend(vspec, rspec):
    from chimera.specifier import evalSpec
    vlist = evalSpec(vspec).models()
    from VolumeViewer import Volume
    vlist = [v for v in vlist if isinstance(v,Volume)]
    res = evalSpec(rspec).models()
    from _volume import Volume_Model
    res = [v for v in res if isinstance(v,Volume_Model)]
    if len(res) == 0:
        res = Volume_Model()
        res.name = 'blend'
        from chimera import openModels
        openModels.add([res], baseId = int(rspec.lstrip('#')))
    else:
        res = res[0]
    blend_maps(vlist, res)

if 'arguments' in globals() and len(arguments) == 2:
    blend(*arguments)
else:
    from chimera import replyobj
    replyobj.status('Need two arguments: runscript blend.py #1-3 #5')