schottky_circles_example_b3d_sv.py

"""
in centers_in v
in radii_in s
in max_level s d=2 n=2
out centers_out v
out radii_out s
"""

# Schottky Circles
# Implementation is based on https://github.com/soma-arc/SchottkyLink 
# Refrences: http://klein.math.okstate.edu/IndrasPearls/
import math

class Circle:
    def __init__(self, center, radius):
        self.center = center
        self.radius = radius

    def invertOnPoint(self, point):
        return (point - self.center) * ((self.radius * self.radius) / \
                (np.linalg.norm(point - self.center) * \
                np.linalg.norm(point - self.center))) + self.center

    def invertOnCircle(self, circle):
        coeffR = circle.radius * math.sqrt(2.0) / 2.0
        p1 = self.invertOnPoint(circle.center + np.array([coeffR, coeffR, 0]))
        p2 = self.invertOnPoint(circle.center + np.array([-coeffR, -coeffR, 0]))
        p3 = self.invertOnPoint(circle.center + np.array([coeffR, -coeffR, 0]))
        return Circle._fromPoints(p1, p2, p3)
        
    @classmethod
    def _fromPoints(cls, a, b, c):
        lA = np.linalg.norm(b - c)
        lB = np.linalg.norm(a - c)
        lC = np.linalg.norm(a - b)
        coeffA = lA * lA * (lB * lB + lC * lC - lA * lA)
        coeffB = lB * lB * (lA * lA + lC * lC - lB * lB)
        coeffC = lC * lC * (lA * lA + lB * lB - lC * lC)
        denom = coeffA + coeffB + coeffC
        center = np.array([(coeffA * a[0] + coeffB * b[0] + coeffC * c[0]) / denom, \
                    (coeffA * a[1] + coeffB * b[1] + coeffC * c[1]) / denom, 0])
        radius = np.linalg.norm(center - a)
        return Circle(center, radius)

centers_out = []
radii_out = []
for centers, radii in zip(centers_in, radii_in):
    source_circles = []
    tag = []
    dest_circles = []
    for i, (c, r) in enumerate(zip(centers, radii)):
        original_circle = Circle(np.array(c), r)
        source_circles.append(original_circle)
        dest_circles.append(original_circle)
        tag.append(i)

    num = [0, 4]
    for level in range(0, max_level):
        for i in range(num[level], num[level+1]):
            for target_index in range(0, len(source_circles)):
                if tag[i] == target_index:
                    continue
                circle = Circle(dest_circles[i].center, dest_circles[i].radius)
                newCircle = source_circles[target_index].invertOnCircle(circle)
                
                tag.append(target_index)
                dest_circles.append(newCircle)
                
        num.append(len(dest_circles))

    centers_out.append([c.center.tolist() for c in dest_circles])
    radii_out.append([float(c.radius) for c in dest_circles])