Examples¶

Here are some simple usage examples of the package. The examples can be found in <tt>examples</tt> directory under main directory.

The line:

sys.path.append('../')

in each example is to allow running script from examples directory without need to setup PATH or PYTHON_PATH system variable to contain dmga main directory.

Using Diagrams¶

This example demonstrates basic usage of Power Diagrams, iteration process, statistics collection and simple rendering using draw submodule.

import sys
sys.path.append('../') # for technical reasons

from pydmga.geometry import OrthogonalGeometry
from pydmga.container import Container
from pydmga.diagram import Diagram
from pydmga.draw import assets
from pydmga.draw import render
from random import random
from random import seed
from pydmga.draw.render import Color

# generate random data
(box_x, box_y, box_z) = (10.0, 10.0, 10.0);
count = 8;    
data = [ (i, box_x * random(), box_y * random(), box_z * random(), random()) for i in range(count)]

# setup geometry - periodic in all dimensions
geometry  = OrthogonalGeometry(10, 10, 10, True, True, True);
# setup container for particles
container = Container(geometry)
# add data
container.add(data)   
# create diagram and make it cache cells (cache_on = True)
diagram = Diagram(container, True)
# create renderer to render 3D graphics
display = render.Render(Color.WHITE_OPAQUE)

# to accumulate volume from cells
volume = 0.0
# to restrict possible colors
colors = [(1,0,0,0.3), (0,1,0,0.3), (0,0,1,0.3), (1,1,0,0.3), (0,1,1,0.3), (1,0,1,0.3)]
# iterate over all cells
for i, cell in enumerate(diagram):
    print "Entering cell", i
    (id, x, y, z, r) = container.get(i)
    # add sphere on atom
    display.add(render.Asset(assets.sphere, (x, y, z), r, Color.BLUE_QUARTER))
    # accumulate volume
    volume += cell.volume()    
    neighbours = []
    all_area = 0
    # iterate over all sides of the Voronoi Cell
    for j, side in enumerate(cell.sides):
        # get the side as a list of coordinates in 3D space
        polygon = side.as_coords()
        # get and remember neighbour of the atom, that created this side
        neighbours.append(side.neighbour)                
        # add outline of the side
        display.add(render.Outline(polygon))
        # add side in a choosen color
        display.add(render.ConvexPolygons(polygon, colors[i % len(colors)]))
        # accumulate side area
        all_area += side.area()
    print "Sum of area sides is", all_area, "should be", cell.area()
    print "Volume of cell is", cell.volume()
    print "Neighbours of cell are", neighbours
    print "Leaving cell", i

# compare (for test) volume of the box and the accumulated volume of the cells
print "Volume is", volume, "should be", box_x * box_y * box_z
# show graphics
display.run()

Using Projective Geometries¶

This example demonstrates basic usage of Projective Geometries. Projective Geometries are Geometry objects that transforms data before insertion to the container. Usually it is some 2D surface of some shape in 3D space. For example, in CastCylinderGeometry the particles are cast onto the side of a given cylinder. Then the 3D Power Diagram of a given shape is computed. Then one can access the geometry of the surface by testing sides of the Voronoi Cells by the geometry on_boundary() function. Function on_boundary() returns True if the cell side is intersection of cell with the 2D surface. Thus we can create the (approximate) Voronoi diagram on 2D surface.

import sys
sys.path.append('../') # for technical reasons

from pydmga.geometry import CastCylinderGeometry
from pydmga.container import Container
from pydmga.diagram import Diagram
from pydmga.draw import assets
from pydmga.draw import render
from pydmga.draw.render import Asset
from pydmga.draw.render import Color
from random import random
from random import seed
from math import pi

# create some random data inside a cylinder
(R, H) = (5.0, 10.0);
count = 64;    
data = [ (i, 2.0 * R * random() - R, 2.0 * R * random() - R, H * random(), 0.1) for i in range(count)]

# create projective geometry
geometry = CastCylinderGeometry(R, H, 128);   
container = Container(geometry)
# add data
container.add(data)
# create Power Diagram
diagram = Diagram(container)
# make renderer
display = render.Render(Color.WHITE_OPAQUE)

colors = [(1,0,0,0.3), (0,1,0,0.3), (0,0,1,0.3), (1,1,0,0.3), (0,1,1,0.3), (1,0,1,0.3)]
all_area = 0
# iterate over all cells
for i, cell in enumerate(diagram):   
    (id, x, y, z, r) = container.get(i)
    display.add(render.Asset(assets.sphere, (x, y, z), r, Color.BLUE_QUARTER)) 
    for j, side in enumerate(cell.sides):
        # if the side lies on the surface
        if (geometry.on_boundary(side.neighbour)):
            polygon = side.as_coords()
            # add polygon twice (once reversed) to allow rendering of both sides
            display.add(render.ConvexPolygons(polygon, colors[i % len(colors)]))
            display.add(render.ConvexPolygons(reversed(polygon), colors[i % len(colors)]))
            # accumulate area
            all_area += side.area()

# adjust area allows to get exact results from inexact computations
print "All area is", geometry.adjust_area(all_area), "should be", (2.0 * pi * R * H)
# show result
display.run()

Using Alpha Shapes¶

Alpha Shape is usually computed using Regular triangulation - a dual to the Power Diagram. In DMG-alpha we assign alpha-values straight to the elements of Power Diagram assuming the duality. Interpretation of the results is almost the same - if alpha-value for some element is smaller than choosen parameter alpha then the element is at least partly unavailable for sphere with radius alpha (in the sense of Laguerre geometry). This allows for detection of cavites and pores by checking for not ocluded vertices and edges in Power Diagram.

import sys
sys.path.append('../') # for technical reasons

from pydmga.geometry import OrthogonalGeometry
from pydmga.container import Container
from pydmga.diagram import Diagram
from pydmga.shape import AlphaShape
from pydmga import draw
from pydmga.draw import assets
from pydmga.draw import Asset
from pydmga.draw import render
from pydmga.draw.render import Color
from random import random
from random import seed
from math import sqrt, pi, sin, cos

# generate random data
(box_x, box_y, box_z) = (10.0, 10.0, 10.0);
count = 10;    
data = [ (i, box_x * random(), box_y * random(), box_z * random(), random()) for i in range(count)]

# setup container for atoms
geometry  = OrthogonalGeometry(10, 10, 10, True, True, True);
container = Container(geometry)
# add data
container.add(data)   

# create basic diagram
diagram = Diagram(container, True)
# create alpha complex
shape = AlphaShape(diagram)

# setup renderer
display = render.Render(Color.WHITE_OPAQUE)
display.show_shape = True

# this is the value at which all elements are in the alpha-shape
# the alpha shape then is topologically equivalent to convex hull of all points
max_alpha = shape.max_alpha_threshold()
# for test
alpha_for_shape = max_alpha / 3.0
print "Max alpha for whole thing is", max_alpha, "chosen alpha is", alpha_for_shape
# here we will add vertices to be drawn
vertices = draw.PointCloud()
# iterate over all shape elements - Voronoi cells enchanced with alpha-shape information
for i, cell_shape in enumerate(shape):
    print "max alpha for cell", i, "is", cell_shape.max_alpha_threshold()
    (id,x,y,z,r) = container.get(i)
    # add the atom sphere to the renderer
    display.add(Asset(assets.sphere, (x,y,z), sqrt(r*r + alpha_for_shape), Color.BLUE_QUARTER))
    # we can use the same functions as in case of normal cell...
    for side in cell_shape:
        print side.as_coords()
    # ...but we also can access alpha-shape description
    # if shape() is called without arguments we get all information
    for simplex in cell_shape.shape():
        # cast information to simpler formats
        print "simplex:"
        print "\t", simplex.as_tuple()
        print "\t", simplex.as_dict()
        print "\t", simplex.as_dict_short()
        # alpha_threshold is an alpha value at which this element belongs to alpha-shape
        # we want to draw only those elements that are greater than choosen alpha_for_shape
        # we draw them thin and in gray color 
        if (simplex.alpha_threshold > alpha_for_shape):
            color = Color.GRAY30_QUARTER
            if   simplex.dimension == 1: 
                part = simplex.diagram_part()
                print "Side is:", part.vertices
            elif simplex.dimension == 2: 
                part = simplex.diagram_part()
                print "Edge is:", part.as_pair()
                display.add(draw.Outline(part.as_coords(), 1.0, color))
            elif simplex.dimension == 3: 
                part = simplex.diagram_part()
                print "Vertex is:", part.id, ",", part.as_coords()
                vertices.add_vertex(part.as_coords(), 2.0, color)
    # next, we iterate over all elements in the concrete shape for given alpha_for_shape parameter
    for simplex in cell_shape.shape(alpha_for_shape):
        print "simplex:"
        print "\t", simplex.as_tuple()
        print "\t", simplex.as_dict()
        print "\t", simplex.as_dict_short()
        # we draw the thing with green color and thick
        color = Color.GREEN_QUARTER
        if   simplex.dimension == 1:
            part = simplex.diagram_part()
            # we add it twice (one time reversed) to get both sides rendered
            display.add(draw.ConvexPolygons(part.as_coords(), color))
            display.add(draw.ConvexPolygons(reversed(part.as_coords()), color))
        elif simplex.dimension == 2: 
            display.add(draw.Outline(simplex.diagram_part().as_coords(), 5.0, color))
        elif simplex.dimension == 3: 
            vertices.add_vertex(simplex.diagram_part().as_coords(), 15.0, color)
# add all vertices at once to the renderer
display.add(vertices)
# run rendering
display.run()

Using SASA Shapes¶

Solvent Accessible Surface is a surface of an union of balls with Van der Waals radii centered at atom positions. Usually it is computed approximately with help of Monte Carlo methods but Power Diagrams may be used to infer the surface in rigorous and exact way. DMG-alpha computes Solvent Accessible Surface Area using spherical geometry - by virtually creating the contour of ocluded sphere surfaces and then computing its area. The difference between sum of all spheres areas and sum of all ocludded areas is an Solvent Accessible Surface Area for a given set of spheres.

import sys
sys.path.append('../') # for technical reasons

from pydmga.geometry import OrthogonalGeometry
from pydmga.container import Container
from pydmga.diagram import Diagram
from pydmga.shape import SASAShape
from pydmga import draw
from pydmga.draw import Asset
from pydmga.draw import assets
from pydmga.draw import render
from pydmga.draw.render import Color
from random import random
from random import seed
from math import pi

# generate some data
(box_x, box_y, box_z) = (10.0, 10.0, 10.0);
count = 5;
data = [(5,5,5,2), (5-2,5,5,2), (5,5-2,5,2), (5,5,5-2,2), (5,5,5+2,2)]

geometry  = OrthogonalGeometry(box_x, box_y, box_z, True, True, True);
container = Container(geometry)
container.add(data)   

# make basic diagram
diagram = Diagram(container)
# create SAS Shape
shape = SASAShape(diagram)

display = render.Render(Color.WHITE_OPAQUE)

# show balls and its periodic images for better understanding
area = shape.sas_area()
print "SASA Area:", area
colors = [Color.RED_OPAQUE]
for ball in container:
    (id, x, y, z, r) = ball
    display.add(Asset(assets.sphere, (x,y,z), r, Color.BLUE_QUARTER))
    display.add(Asset(assets.sphere, (x+box_x,y,z), r, (1.0, 1.0, 1.0, 0.1)))
    display.add(Asset(assets.sphere, (x-box_x,y,z), r, (1.0, 1.0, 1.0, 0.1)))
    display.add(Asset(assets.sphere, (x,y+box_y,z), r, (1.0, 1.0, 1.0, 0.1)))
    display.add(Asset(assets.sphere, (x,y-box_y,z), r, (1.0, 1.0, 1.0, 0.1)))
    display.add(Asset(assets.sphere, (x,y,z+box_z), r, (1.0, 1.0, 1.0, 0.1)))
    display.add(Asset(assets.sphere, (x,y,z-box_z), r, (1.0, 1.0, 1.0, 0.1)))

# get the SASA shape of the first atom
i = 0
cell = shape.get_cell(i)
# render sphere
(id, x, y, z, r) = container.get(i)
display.add(Asset(assets.sphere, (x,y,z), r, Color.BLUE_QUARTER))
# get all area (should be equal to the ball area)
all_area = cell_shape.sas_area() + cell_shape.excluded_area()    
print "area is", cell_shape.sas_area(), ", excluded area is", cell_shape.excluded_area()
print "sum is", all_area, "should be", 4.0 * pi * (r**2);
# you can also access standard functions for basic Voronoi cell
print cell_shape.area(), ",", cell_shape.volume()
# render sides of neighbours other than 0 (0 may be neighbour of 0 since we have Periodic Boundary Conditions)
b = 0
for a, side in enumerate(cell_shape):
    if side.neighbour >= 0:
        b += 1
        display.add(render.Outline(side.as_coords(), 2.0, Color.BLACK_OPAQUE))
# draw outline of the excluded surface
points = draw.PointCloud()
great_arcs = render.GreatArcs((5,5,5),1.0,2.0,Color.BLUE_OPAQUE)
for j, arc in enumerate(cell_shape.shape()):       
    print "Arc is: ", arc
    arcs = draw.SphereArcs((x,y,z), r, 2.0, colors[j % len(colors)])
    arcs.add_arc(arc.on_plane, arc.first, arc.second)
    great_arcs.add_arcs([arc.first, arc.on_sphere, arc.second], False)
    points.add_vertex(arc.first, 10.0, Color.GREEN_OPAQUE)
    points.add_vertex(arc.second, 10.0, Color.GREEN_OPAQUE)
    display.add(arcs)
display.add(great_arcs)
display.add(points)

display.run()

Drawing¶

DMG-alpha uses Panda3D rendering engine to provide simple drawing in 3D. panda3D can be downloaded from http://panda3d.org .

import sys
sys.path.append('../') # for technical reasons

from pydmga import draw
from pydmga.draw import render
from pydmga.draw import assets
from math import sqrt
import random

# create renderer object
r = draw.Render(render.Color.WHITE_OPAQUE)
# add convex polygon in 3D (you must actually pass something that is covex and is polygon, otherwise bad things may happen)
r.add(draw.ConvexPolygons([(-1.0, -1.0, 0.0),(-1.0, 1.0, 0.0),(1.0, 1.0, 0.5)], (1.0, 0.0, 0.0, 0.5)))
r.add(draw.ConvexPolygons([(-1.0, -1.0, 0.0),(1.0, 1.0, 0.5),(-1.0, 1.0, 0.0)], (1.0, 0.0, 0.0, 0.5)))
# add collecion of lines (a triangle in this example) with line width equal to 2.0 and a given color. 
# by default Outline creates closed contour but you can turn this off with additional parameter set to False
r.add(draw.Outline([(-1.0, -1.0, 0.0),(1.0, 1.0, 0.5),(-1.0, 1.0, 0.0)], 2.0, (1.0, 1.0, 0.0, 0.5)))

# make some random points and place them on the sphere
v = [(random.random() * 10.0 - 5.0, random.random() * 10.0 - 5.0, random.random() * 10.0 - 5.0) for i in range(100)]
for i in range(len(v)):
    d = sqrt((v[i][0])**2 + (v[i][1])**2 + (v[i][2])**2);
    if (d > 5.0):
        v[i] = (5.0 * v[i][0] / d, 5.0 * v[i][1] / d, 5.0 * v[i][2] / d)

# add points to renderer, with point size equal 5.0 and with given color
r.add(draw.PointCloud().add_vertices(v, 5.0, (0.5,0.5,0.5,1.0)))
# add a sphere with BLUE color and 75% transparent (25% opaque)
# you can use other predefined colors: RED, BLUE, GREEN, GRAY, WHITE, BLACK
# you can use weaker and stronger colors by adding number 10, 20, .., 90 to the name of the color
# e.g. GREEN50 is darker color than GREEN90. GREEN is the name for full GREEN i.e. rgb = (0, 255, 0)
# each color has three variants: OPAQUE, HALF, QUARTER (opaque, half-transparent, 75% transparent)
r.add(draw.Asset(assets.sphere, (0.0, 0.0, 0.0), 5.0, render.Color.BLUE_QUARTER))
# draws great arcs on the side of the sphere
# constructor takes parameters of the sphere
# then add_arcs() add consecutive arcs from one point to another
# by default the shape is closed, but this can be turned off
ga = draw.GreatArcs((0.0, 0.0, 0.0), 5.0).add_arcs([(5.0,0.0,0.0),(0.0,5.0,0.0),(0.0,0.0,5.0)])
r.add(ga)
# draws points on the joins of the arcs
r.add(draw.PointCloud().add_vertices(ga.vertices, 12.0, (0.0, 1.0,0.0,1.0)))
# creates cutting plane arc - an intersection of sphere and some polygon
r.add(draw.SphereArcs((0, 0, 0), 5, 2, (0,1,0,1)).add_cutting_plane_arc((-1, 0.0, 0.0), 4.0, (-4,0,3),(-4,3,0)))
r.add(draw.SphereArcs((0, 0, 0), 5, 2, (0,0,1,1)).add_cutting_plane_arc((-1, 0.0, 0.0), 4.0, (-4,3,0),(-4,0,3)))
r.add(draw.SphereArcs((0, 0, 0), 5, 2, (0,1,0,1)).add_cutting_plane_arc((0, 1.0, 0.0), 4.0, (0,4,3),(0,4,-3)))
r.add(draw.SphereArcs((0, 0, 0), 5, 2, (0,0,1,1)).add_cutting_plane_arc((0, 1.0, 0.0), 4.0, (0,4,-3),(0,4,3)))
r.add(draw.SphereArcs((0, 0, 0), 5, 2, (1,0,0,1)).add_cutting_plane_arc((0, 0.0, 1.0), 0.0, (0,-5,0),(-5,0,0)))
r.add(draw.SphereArcs((0, 0, 0), 5, 2, (0.5,1,0,1)).add_cutting_plane_arc((0, 0.0, 1.0), 0.0, (-5,0,0), (0,-5,0)))
# add vertices on the joint points
r.add(draw.PointCloud().add_vertices([(-4,0,3),(-4,3,0),(0,4,-3),(0,4,3)], 8.0, (0.0, 1.0,1.0,1.0)))
r.add(draw.PointCloud().add_vertices([(0,-5,0),(-5,0,0)], 8.0, (1.0, 0.0,0.0,1.0)))
# add some cube
r.add(draw.Asset(assets.cube, (-5.0, 5.0, 5.0), 1.0, render.Color.GREEN_OPAQUE))
# run rendering
r.run()