"""Define algorithms to support intersection calculation."""
__author__ = "Daniel Ching, Doga Gursoy"
__copyright__ = "Copyright (c) 2016, UChicago Argonne, LLC."
__docformat__ = 'restructuredtext en'
__all__ = [
'clip_SH',
'halfspacecirc',
]
import logging
from math import asin, sqrt
import warnings
import numpy as np
from xdesign.geometry.point import *
logger = logging.getLogger(__name__)
[docs]def halfspacecirc(d, r):
"""Return the area of intersection between a circle and half-plane.
Returns the smaller fraction of a circle split by a line d units away
from the center of the circle.
Parameters
----------
d : scalar
The distance from the line to the center of the circle
r : scalar
The radius of the circle
Returns
-------
f : scalar
The proportion of the circle in the smaller half-space
References
----------
Glassner, A. S. (Ed.). (2013). Graphics gems. Elsevier.
"""
assert r > 0, "The radius must positive"
assert d >= 0, "The distance must be positive or zero."
if d >= r: # The line is too far away to overlap!
return 0
f = 0.5 - d * sqrt(r**2 - d**2) / (np.pi * r**2) - asin(d / r) / np.pi
# Returns the smaller fraction of the circle, so it can be at most 1/2.
if f < 0 or 0.5 < f:
if f < -1e-16: # f will often be less than 0 due to rounding errors
warnings.warn(
"halfspacecirc was out of bounds, {}".format(f), RuntimeWarning
)
f = 0
return f
def two_lines_intersect(l0A, l0b, l1A, l1b):
A = np.stack([l0A, l1A], axis=0)
b = np.stack([l0b, l1b], axis=0)
x = np.linalg.solve(A, b)
return x
def half_space(self, center):
"""Return the half space polytope respresentation of the Line."""
A, B = self.standard
# test for positive or negative side of line
if not halfspace_has_point(A, B, center):
A = -A
B = -B
return A, B
def halfspace_has_point(A, B, point):
return np.dot(A, point._x) <= B
[docs]def clip_SH(clipEdges, polygon):
"""Clip a polygon using the Sutherland-Hodgeman algorithm.
Parameters
----------
clipEdges [[A, b], ...]
half-spaces defined by coefficients
polygon
"""
outputList = polygon.vertices
for clipEdge in clipEdges:
# previous iteration output is this iteration input
inputList = outputList
outputList = list()
if len(inputList) == 0:
break
S = inputList[-1]
for E in inputList:
if halfspace_has_point(clipEdge[0], clipEdge[1], E):
if not halfspace_has_point(clipEdge[0], clipEdge[1], S):
A, b = calc_standard(np.stack([S._x, E._x], axis=0))
new_vert = two_lines_intersect(
A, b, clipEdge[0], clipEdge[1]
)
outputList.append(Point(new_vert))
outputList.append(E)
elif halfspace_has_point(clipEdge[0], clipEdge[1], S):
A, b = calc_standard(np.stack([S._x, E._x], axis=0))
new_vert = two_lines_intersect(A, b, clipEdge[0], clipEdge[1])
outputList.append(Point(new_vert))
S = E
return outputList