# -*- coding: utf-8 -*-
"""
Created on Mon May 25 15:23:49 2015
@author: chris
based on add_isosurface function from chemview
"""
import numpy as np
#TODO not making this a depedency until it works
try:
from skimage.measure import marching_cubes, correct_mesh_orientation
except ImportError:
pass
[docs]def get_isosurface(coordinates, function, isolevel=0.03, color=(255, 0, 0, 255),
extents=(5,5,5), resolution=100):
'''Add an isosurface to the current scene.
Parameters
----------
coordinates : numpy.array
the coordinates of the system
function : func
A function that takes x, y, z coordinates as input and is broadcastable
using numpy. Typically simple functions that involve standard
arithmetic operations and functions such as ``x**2 + y**2 + z**2`` or
``np.exp(x**2 + y**2 + z**2)`` will work. If not sure, you can first
pass the function through ``numpy.vectorize``.
Example: ``mv.add_isosurface(np.vectorize(f))``
isolevel : float
The value for which the function should be constant.
color : (int, int, int, int)
The color given in RGBA format
extents : (float, float, float)
+/- x,y,z to extend the molecule geometrty when constructing the surface
resolution : int
The number of grid point to use for the surface. An high value will
give better quality but lower performance.
'''
if coordinates.shape[0] == 0:
return False
# We want to make a container that contains the whole molecule
# and surface
coordinates = coordinates.astype('float32')
area_min = coordinates.min(axis=0) - np.array(extents)
area_max = coordinates.max(axis=0) + np.array(extents)
x = np.linspace(area_min[0], area_max[0], resolution)
y = np.linspace(area_min[1], area_max[1], resolution)
z = np.linspace(area_min[2], area_max[2], resolution)
xv, yv, zv = np.meshgrid(x, y, z)
spacing = np.array((area_max - area_min)/resolution)
if isolevel >= 0:
sign = 1
else:
sign = -1
volume = sign*function(xv, yv, zv)
if volume.flatten().min() <= sign*isolevel and volume.flatten().max() >= sign*isolevel:
verts, faces = marching_cubes(volume, sign*isolevel)
#don't know if i need this
faces = correct_mesh_orientation(volume, verts, faces,
gradient_direction='descent')
verts = area_min + spacing/2 + verts*spacing
# TODO for some reason need to invert, but no one knows what the problem is
verts[:,[0,1]] = verts[:,[1,0]]
# TODO and its messing up the normals so using double facing
#(kind of works if transparent)
opposite_faces = faces.copy()
opposite_faces[:,[0,1]] = opposite_faces[:,[1,0]]
faces = np.concatenate([faces,opposite_faces])
else:
return
normals = calc_normals(verts, faces)
verts = verts[faces.flatten()]
normals = normals[faces.flatten()]
colors = np.array([color for _ in range(verts.shape[0])])
return verts.astype('float32'), normals.astype('float32'), colors
[docs]def normalize_v3(arr):
''' Normalize a numpy array of 3 component vectors shape=(n,3) '''
lens = np.sqrt( arr[:,0]**2 + arr[:,1]**2 + arr[:,2]**2 )
arr[:,0] /= lens
arr[:,1] /= lens
arr[:,2] /= lens
return arr
[docs]def calc_normals(verts, faces):
"""from;
https://sites.google.com/site/dlampetest/python/calculating-normals-of-a-triangle-mesh-using-numpy
"""
#Create a zeroed array with the same type and shape as our vertices i.e., per vertex normal
norm = np.zeros( verts.shape, dtype=verts.dtype )
#Create an indexed view into the vertex array using the array of three indices for triangles
tris = verts[faces]
#Calculate the normal for all the triangles, by taking the cross product of the vectors v1-v0, and v2-v0 in each triangle
n = np.cross( tris[::,1 ] - tris[::,0] , tris[::,2 ] - tris[::,0] )
# n is now an array of normals per triangle. The length of each normal is dependent the vertices,
# we need to normalize these, so that our next step weights each normal equally.
normalize_v3(n)
# now we have a normalized array of normals, one per triangle, i.e., per triangle normals.
# But instead of one per triangle (i.e., flat shading), we add to each vertex in that triangle,
# the triangles' normal. Multiple triangles would then contribute to every vertex, so we need to normalize again afterwards.
# The cool part, we can actually add the normals through an indexed view of our (zeroed) per vertex normal array
norm[ faces[:,0] ] += n
norm[ faces[:,1] ] += n
norm[ faces[:,2] ] += n
normalize_v3(norm)
return norm
if __name__ == '__main__':
from chemlab.io import remotefile
url = "https://raw.githubusercontent.com/cclib/cclib/master/data/GAMESS/basicGAMESS-US2012/water_mp2.out"
df = remotefile(url, "gamess")
molecule = df.read('molecule')
molecule.guess_bonds()
mocoeffs = df.read('mocoeffs')
gbasis = df.read('gbasis')
from chemlab.qc import molecular_orbital
orbital = 2
iso_level = -0.3
color = (255, 0, 0, 100)
coefficients = mocoeffs[0][orbital]
f = molecular_orbital(molecule.r_array, coefficients, gbasis)
verts, normals, colors = get_isosurface(molecule.r_array, f,
iso_level, color)