Skip to content

Instantly share code, notes, and snippets.

@jasonmhite
Last active January 12, 2018 19:06
Show Gist options
  • Save jasonmhite/b61fdac3ec38df9323a329d73ff602f9 to your computer and use it in GitHub Desktop.
Save jasonmhite/b61fdac3ec38df9323a329d73ff602f9 to your computer and use it in GitHub Desktop.
Computation of the solid angle of an oriented right triangle in free space
##################################################################################
# Copyright (c) 2018, Jason M. Hite
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# 1. Redistributions of source code must retain the above copyright notice, this
# list of conditions and the following disclaimer.
# 2. Redistributions in binary form must reproduce the above copyright notice,
# this list of conditions and the following disclaimer in the documentation
# and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
# ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
##################################################################################
# This file gives an example of how to calculate the solid angle of an arbitrarily
# oriented right-triangle as viewed from a specified vantage point. It implements
# the method of Van Oosteram and Strackee in an efficient vectorized form using
# NumPy.
#
# Refs
# Van Oosterom, A., & Strackee, J. (1983). The solid angle of a plane triangle.
# IEEE transactions on Biomedical Engineering, (2), 125-126.
#
# A fast scalar triple product in NumPy
# https://stackoverflow.com/questions/20908754/how-to-speed-up-a-vector-cross-product-calculation/20910319#20910319
#
# Basic Guide to einsum
# http://ajcr.net/Basic-guide-to-einsum/
import numpy as np
# Levi-Civita symbol
eijk = np.zeros((3, 3, 3))
eijk[0, 1, 2] = eijk[1, 2, 0] = eijk[2, 0, 1] = 1
eijk[0, 2, 1] = eijk[2, 1, 0] = eijk[1, 0, 2] = -1
def triangle_sa(R1, R2, R3):
"""
Compute the solid angle of an arbitrary triangle. R1, R2, R3 are
3-vectors from the vantage point to the vertices of the triangle.
I.e., w.l.o.g. the vantage point is assumed to lie at the origin.
"""
l = np.linalg.norm(np.vstack((R1, R2, R3)), axis=1)
p = np.array([R2.dot(R3), R1.dot(R3), R1.dot(R2)])
N = np.einsum('ijk,i,j,k->', eijk, R1, R2, R3) # scalar triple
D = l.prod() + l.dot(p)
return 2. * np.arctan2(N, D)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment