triangle_loci.py

from PIL import Image, ImageDraw, ImageSequence
from random import random as rnd
from math import sqrt

def vec3_scale(v,f):
    return (v[0]*f,v[1]*f,v[2]*f)

def vec3_add(a,b):
    return (a[0]+b[0],a[1]+b[1],a[2]+b[2])

def vec3_sub(a,b):
    return (a[0]-b[0],a[1]-b[1],a[2]-b[2])

def vec3_dot(a,b):
    return a[0]*b[0]+a[1]*b[1]+a[2]*b[2]

def vec3_cross(a,b):
    return (a[1]*b[2]-a[2]*b[1],a[2]*b[0]-a[0]*b[2],a[0]*b[1]-a[1]*b[0])
    
def vec3_mag(v):
    return sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2])
    
#// Compute barycentric coordinates (u, v, w) for
#// point p with respect to triangle (a, b, c)
#void Barycentric(Point a, Point b, Point c, float &u, float &v, float &w)
#{
#    Vector v0 = b - a, v1 = c - a, v2 = p - a;
#    float d00 = Dot(v0, v0);
#    float d01 = Dot(v0, v1);
#    float d11 = Dot(v1, v1);
#    float d20 = Dot(v2, v0);
#    float d21 = Dot(v2, v1);
#    float denom = d00 * d11 - d01 * d01;
#    v = (d11 * d20 - d01 * d21) / denom;
#    w = (d00 * d21 - d01 * d20) / denom;
#    u = 1.0f - v - w;
#}


BLACK = (0,0,0)
GREEN = (0,255,0)
BLUE = (0,0,255)
RED = (255,0,0)
WHITE = (255,255,255)
CYAN = (128,255,255)
YELLOW = (255,255,128)

def tri_ray_loci(a,b,c,ray_origin,ray_dir,ray_radius):
    # get triangle edge vectors and plane normal
    u = vec3_sub(b,a)
    v = vec3_sub(c,a)
    n = vec3_cross(u,v)
    if n[0]==0 and n[1]==0 and n[2]==0:
        return RED # triangle is degenerate
    w0 = vec3_sub(ray_origin,a)
    j = vec3_dot(n,ray_dir)
    if abs(j) < 0.00000001:
        #TODO parallel?
        return BLACK # parallel, disjoint or on plane
    # get intersect point of ray with triangle plane
    i = -vec3_dot(n,w0)
    k = i / j
    #if k < 0.0:
    #    return AWAY # ray goes away from triangle
    # for a segment, also test if (k > 1) => no intersect
    hit = vec3_add(ray_origin,vec3_scale(ray_dir,k)) # intersect point of ray and plane
    # barycentric coords
    w = vec3_sub(hit,a)
    uu = vec3_dot(u,u)
    uv = vec3_dot(u,v)
    uw = vec3_dot(u,w)
    vv = vec3_dot(v,v)
    vw = vec3_dot(v,w)
    invDenom = 1.0 / (uu * vv - uv * uv)
    U = (vv * uw - uv * vw) * invDenom
    V = (uu * vw - uv * uw) * invDenom
    W = U+V
    # basic in-triangle test
    A, B, C = U>=0.0, V>=0.0, W<1.0
    if A and B and C:
        return GREEN
    # near miss?
    if U>-0.1 and B and C:
        return BLUE
    if A and V>-0.1 and C:
        return YELLOW
    if A and B and W<1.1:
        return CYAN
    # far miss
    return BLACK

def test():
    dim = 500
    a = (rnd()*dim,rnd()*dim,rnd()*dim)
    b = (rnd()*dim,rnd()*dim,rnd()*dim)
    c = (rnd()*dim,rnd()*dim,rnd()*dim)
    img = Image.new("RGB",(dim,dim))
    pix = img.load()
    for y in xrange(dim):
        for x in xrange(dim):
            col = tri_ray_loci(a,b,c,[x,y,0],[0,0,dim],30)
            pix[x,y] = col
    def marker(x,y,col):
        for i in range(-1,2):
            for j in range(-1,2):
                pix[round(x)+j,round(y)+i] = col
    marker(a[0],a[1],WHITE)
    marker(b[0],b[1],WHITE)
    marker(c[0],c[1],WHITE)
    img.show()
    
test()