geo diff

gistfile1.py

import numpy
import pprint
import random as rndom
from plyfile import PlyData# PlyElement
from scipy import spatial  # kd tree
from collections import defaultdict

######UTILS#######
def get_coordinates(vertex_id, mesh):
    vertex = mesh['vertex'].data[vertex_id]
    return numpy.array((vertex['x'], vertex['y'], vertex['z']))

def vec_ang(v1, v2):
    """ Returns the angle in radians between vectors 'v1' and 'v2'  based on https://newtonexcelbach.wordpress.com/2014/03/01/the-angle-between-two-vectors-python-version/  """
    cosang = numpy.dot(v1, v2)
    sinang = numpy.linalg.norm(numpy.cross(v1, v2))
    return numpy.arctan2(sinang, cosang)

def check_1_ring(neighbours):
    """ just a sanity test to make sure i've implemented 1-ring correctly """
    for key in neighbours:
        for val in neighbours[key]:
            if key not in neighbours[val]:
                print val, " doesnt have ", key, " but ", key, " does "


def generate_1_ring(mesh):
    """generates a dictionary mapping each vertext to a list of its 1-ring neighbours
    TODO: why is dict len less than number of vertices?"""
    neighbours = defaultdict(set)
    for face_indices in mesh['face'].data['vertex_index']:
        for v_index in face_indices:
            f_ind = face_indices.tolist()
            f_ind.remove(v_index)
            neighbours[v_index].update(f_ind)
    return dict(neighbours)

def generate_face_mapping(mesh):
    """generate the mapping from a vertice to the faces it is a member of"""
    mapping = defaultdict(list)
    for i,face_indices in enumerate(mesh['face'].data['vertex_index']):
        for v_index in face_indices:
            mapping[v_index].append(i)
    return dict(mapping)
######### END UTILS #######


def uniform_laplace(vertex_id, neighbours, mesh):
    if vertex_id not in neighbours:
        raise Exception("it hath no neighbours")

    summation = numpy.array([0, 0, 0])
    vertex = get_coordinates(vertex_id, mesh)
    weighting = 1. / len(neighbours[vertex_id])
    for neighbour in neighbours[vertex_id]:
        summation = summation + (get_coordinates(neighbour, mesh) - vertex)
    return weighting * summation

def generate_uni_laplace(neighbours, mesh):
    uni_laplace = dict()
    for vertex_id in neighbours:
        uni_laplace[vertex_id] = uniform_laplace(vertex_id, neighbours, mesh)
    return uni_laplace

def mean_curvature_from_coords(diff_coord):
    """.5 times absolute value"""
    return .5 * numpy.linalg.norm(diff_coord)

def generate_mean_curvature(uni_laplace):
    """Should this all be bary centered???"""

    mean_curvature = dict()
    for vertex_id in uni_laplace:
        mean_curvature[vertex_id] = mean_curvature_from_coords(uni_laplace[vertex_id])
    return mean_curvature


def angle(a,b,c):
    """returns angle formed between ab (b - a) and ac (c-a)"""
    return vec_ang(b-a, c-a)
def angle_deficit(vertex_c):
    """ computes the angle deficit between vertex and all its 1 ring neighbours"""


def main(mesh=False):
    if mesh is False:  # q is being registered to reference p
        mesh = PlyData.read(open('/Users/md/Documents/UCL 3/geometry processing/cwk2/code/sphere.ply'))
    neighbours = generate_1_ring(mesh)
    lap = generate_uni_laplace(neighbours, mesh)
    mean = generate_mean_curvature(lap)
    pprint.pprint(mean)

if __name__ == "__main__":
    main()