Hmm. The lack of comments make it slightly non-obvious what is happening. The inverted conditions not iszero(…)
don't make it any easier to understand. The same set of these conditions occurs two times, which is a bit confusing. I also don't know where the iszero(…)
function is from, and will replace it by … == 0
in the following.
Now a bit of math first. Two vectors are perpendicular if their scalar product is zero:
1st vector (x, y, z)
2nd vector (a, b, c)
0 = a·x + b·y + c·z
You are correct in avoiding the trivial solution x = y = z = 0
. Note that you do not avoid the solution a = b = c = 0
, because random.random()
can return zero.
If one of x, y, z
is zero, then the above equation can be solved by choosing a non-zero value for the corresponding direction, and setting the other variables to zero. Example: given (1, 0, 3)
as a first vector, the equation can be solved by a second vector (0, 1, 0)
.
The above rule can also be applied in reverse if the first vector has two zero fields, but it's only required to change one of the zero fields to a non-zero value.
We can encode these simple cases in pseudocode as
match vec1 with
| (0, 0, 0) -> ValueError('zero-vector')
| (0, _, _) -> Vector(1, 0, 0)
| (_, 0, _) -> Vector(0, 1, 0)
| (_, _, 0) -> Vector(0, 0, 1)
| (x, y, z) -> a more complex case which we'll handle in a moment
If all parts of the input vector are non-zero, the calculation is a bit more complex because we have three variables a, b, c
to determine but only one equation – this means we can choose two values arbitrarily. Choosing a random value might make sense when hardening your application, but it's difficult to test, and we can make the code simpler by choosing a = b = 1
(we cannot choose a = b = 0
). The equation can now be used to calculate c
:
c = -(x + y)/z
which you essentially have used as well.
So to wrap it all up, I would write this code:
def perpendicular_vector(v):
r""" Finds an arbitrary perpendicular vector to *v*."""
# for two vectors (x, y, z) and (a, b, c) to be perpendicular,
# the following equation has to be fulfilled
# 0 = ax + by + cz
# x = y = z = 0 is not an acceptable solution
if v.x == v.y == v.z == 0:
raise ValueError('zero-vector')
# If one dimension is zero, this can be solved by setting that to
# non-zero and the others to zero. Example: (4, 2, 0) lies in the
# x-y-Plane, so (0, 0, 1) is orthogonal to the plane.
if v.x == 0:
return Vector(1, 0, 0)
if v.y == 0:
return Vector(0, 1, 0)
if v.z == 0:
return Vector(0, 0, 1)
# arbitrarily set a = b = 1
# then the equation simplifies to
# c = -(x + y)/z
return Vector(1, 1, -1.0 * (v.x + v.y) / v.z)