HiveBrain v1.2.0
Get Started
← Back to all entries
patternpythonMinor

Calculate angle between planes

Submitted by: @import:stackexchange-codereview··
0
Viewed 0 times
anglebetweencalculateplanes

Problem

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 +

Solution

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

About precomputing normals:

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

        def calculate_angle():
            ....
            if (....):
                ....
                # Already computed!
                # n1 = angle_norm(index1)
                # n2 = angle_norm(index2)
                ang = angle_between(normals[index1], normals[index2])
                angles_list.append(ang)

Code Snippets

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

        def calculate_angle():
            ....
            if (....):
                ....
                # Already computed!
                # n1 = angle_norm(index1)
                # n2 = angle_norm(index2)
                ang = angle_between(normals[index1], normals[index2])
                angles_list.append(ang)

Context

StackExchange Code Review Q#162624, answer score: 4

Revisions (0)

No revisions yet.