# Calculate angle between planes

I have written working code calculating the angle between the adjacent planes. I read subsequently from standart input: n - amount of triangles, m - amount of vertices ind[] - indices of the vertices given coord[] - coordinates of the vertices given

The output is supposed to be the maximum angle between the adjacent planes in rad.

Function calculate_angle() iterates over the amount of triangles so that I compare only new ones with the old ones (for z in range(k, n))

list_result = [i for i in index1 if i in index2] used for the adjacency detection: it asks, whether the indices of the first triangle coordinates == to the second. If there are at leat 2 of them - we start to calculate the normals to triangles (the angle between surfaces = to the angle between their normals)

However, if the number of triangles is more than 104 it starts to work very slowly:

import numpy as np
import math
import sys
import cProfile, pstats, io
pr = cProfile.Profile()
pr.enable()

num_triangles, num_vertices = input().split()
n = int(num_triangles)  # amount of triangles
m = int(num_vertices)  # amount of vertices
ind = []
coord = []
angles_list =[]

for i in range(n):
ind.append([int(j) for j in input().split()])  # indices of the vertices given
for j in range(m):
coord.append([float(k) for k in input().split()])  # coordinates of the vertices given

def unit_vector(v):
# Returns the unit vector of the vector.
return v/ np.sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2])

def angle_between(v1, v2):
v1_u = unit_vector(v1)
v2_u = unit_vector(v2)
return math.acos(max(min(np.dot(v1_u, v2_u), 1), -1))   # (np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))

def calculate_angle():
for k in range(0, n):
for z in range(k, n):
index1 = ind[k]
index2 = ind[z]
trignum =0
list_result = [i for i in index1 if i in index2]
if (ind[k] != ind[z])&(len(list_result) >= 2)&(trignum <= 3):
trignum = trignum + 1
n1 = angle_norm(index1)
n2 = angle_norm(index2)
ang = angle_between(n1, n2)
angles_list.append(ang)
return max(angles_list)

def cross(a, b):
c = [a[1]*b[2] - a[2]*b[1],
a[2]*b[0] - a[0]*b[2],
a[0]*b[1] - a[1]*b[0]]
return c

def angle_norm(triangle_index):
p0 = coord[triangle_index[0]]
p1 = coord[triangle_index[1]]
p2 = coord[triangle_index[2]]
V1 = np.array(p1) - np.array(p0)
V2 = np.array(p2) - np.array(p0)
an = cross(V1, V2)
return an

print(calculate_angle())

pr.disable()
s = io.StringIO()
sortby = 'tottime'
ps = pstats.Stats(pr, stream=s).sort_stats(sortby)
ps.print_stats()
print(s.getvalue())


The analysis results:

0.0
12130580 function calls in 26.836 seconds

Ordered by: internal time

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
1048576    6.038    0.000    6.038    0.000 /Users/mac/PycharmProjects/untitled3/Course:21(unit_vector)
1284    5.309    0.004    5.309    0.004 {built-in method builtins.input}
4194304    4.645    0.000    4.645    0.000 {built-in method numpy.core.multiarray.array}
1048576    3.304    0.000   10.566    0.000 /Users/mac/PycharmProjects/untitled3/Course:55(angle_norm)
1048576    2.617    0.000    2.617    0.000 /Users/mac/PycharmProjects/untitled3/Course:49(cross)
1    2.090    2.090   21.516   21.516 /Users/mac/PycharmProjects/untitled3/Course:33(calculate_angle)
524288    0.996    0.000    8.224    0.000 /Users/mac/PycharmProjects/untitled3/Course:26(angle_between)
524288    0.531    0.000    0.531    0.000 {built-in method numpy.core.multiarray.dot}
819840    0.487    0.000    0.487    0.000 /Users/mac/PycharmProjects/untitled3/Course:39(<listcomp>)
524288    0.329    0.000    0.329    0.000 {built-in method builtins.min}
524289    0.249    0.000    0.249    0.000 {built-in method builtins.max}
524288    0.095    0.000    0.095    0.000 {built-in method math.acos}
819840    0.080    0.000    0.080    0.000 {built-in method builtins.len}
525571    0.055    0.000    0.055    0.000 {method 'append' of 'list' objects}
1    0.008    0.008    0.008    0.008 {built-in method builtins.print}
1280    0.002    0.000    0.002    0.000 /Users/mac/PycharmProjects/untitled3/Course:16(<listcomp>)
1284    0.000    0.000    0.000    0.000 {method 'split' of 'str' objects}
3    0.000    0.000    0.000    0.000 /Users/mac/PycharmProjects/untitled3/Course:18(<listcomp>)
1    0.000    0.000    0.000    0.000 /Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/codecs.py:318(decode)
1    0.000    0.000    0.000    0.000 {built-in method _codecs.utf_8_decode}
1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
Here's what I already tried to optimise:

1. I got rid of couple of np built-in functions, e.g. np.cross() and np.linalg.norm(), that gave me a couple of seconds.
2. It was for z in range(1, n), I changed 1 to k in order to not take into account already calculated triangles.

I also tried to make faster input, but to no avail. (Used map, used sys.std)

• "I also tried to make faster input, but to no avail." What do you mean in particular, could you extend a bit about that. Also did you already use a decent profiling tool to discover the bottlenecks in your code? – πάντα ῥεῖ May 5 '17 at 17:59
• @πάνταῥεῖ The code clearly uses cProfile. – 200_success May 5 '17 at 18:04
• My bad. I should have pointed that out before directing him here. – Carcigenicate May 5 '17 at 18:04
• It would be nice if you included a small sample of the input in the question. – 200_success May 5 '17 at 18:07
• @πάνταῥεῖ, I added the results into the question. The bottleneck is the calculate_angle() function. I have no better idea than to think i's because of the 2 for-cycles. I think the only faster algorithm would be some kind of a graph -thing, but I'm not an expert, as I pointed out, so I need some help. – TheDoctor May 5 '17 at 18:22

One immediate advice is to precompute angle_norm for each triangle. That alone will give you some boost.

It also lets reformulate the problem as finding a point set diameter. It is a classical problem, and a generic case is quite hard. You may want to look here for an introduction and references.

However, in your case you deal with normalized vectors; their endpoints lay on the unitary sphere, and the problem reduces to a simpler 2D case, which admits the $O(n\log{n})$ solution.

PS

        normals = []
def compute_normals(triangles):
for t in triangles:
normals.append(unit_vector(angle_norm(t))

def calculate_angle():
....
if (....):
....