schottky_spheres_example_b3d_sv.py

"""
in base_center v
in base_radius s
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 Spheres
# Implementation is based on https://github.com/soma-arc/SchottkyLink 
# Refrences: http://klein.math.okstate.edu/IndrasPearls/
import math

class Sphere:
    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

    @classmethod
    def _pivoting(cls, mat, n, k):
        col = k
        maxVal = math.fabs(mat[k][k])
        for i in range(k+1, n):
            if math.fabs(mat[i][k]) > maxVal:
                col = i
                maxVal = math.fabs(mat[i][k])
        if k != col:
            tmp = mat[col]
            mat[col] = mat[k]
            mat[k] = tmp
        return mat
    
    def __makeSphereFromPoints(self, p1, p2, p3, p4):
        p = np.array([p1, p2, p3, p4])

        coefficient = []
        for i in range(0, 3):
            coefficient.append([])
            coefficient[-1].append(2 * (p[i+1][0] - p[i][0]))
            coefficient[-1].append(2 * (p[i+1][1] - p[i][1]))
            coefficient[-1].append(2 * (p[i+1][2] - p[i][2]))
            coefficient[-1].append(-(math.pow(p[i][0], 2) + \
                math.pow(p[i][1], 2) + math.pow(p[i][2], 2)) + \
                math.pow(p[i+1][0], 2) + \
                math.pow(p[i+1][1], 2) + math.pow(p[i+1][2], 2))
                
        n = 3
        for k in range(0, n-1):
            coefficient = Sphere._pivoting(coefficient, n, k)
            
            vkk = coefficient[k][k]
            for i in range(k+1, n):
                vik = coefficient[i][k]
                for j in range(k, n+1):
                    coefficient[i][j] = coefficient[i][j] - \
                    vik * (coefficient[k][j] / vkk)
        
        coefficient[n-1][n] = coefficient[n-1][n] / coefficient[n-1][n-1]
        for i in range(n-2, -1, -1):
            acc = 0.0
            for j in range(i+1, n):
                acc += coefficient[i][j] * coefficient[j][n]
            coefficient[i][n] = (coefficient[i][n] - acc) / \
                coefficient[i][i]
        
        center = np.array([coefficient[0][3], coefficient[1][3], \
                            coefficient[2][3]])
        radius = np.linalg.norm(center - p1)
        return Sphere(center, radius)
        
    def invertOnSphere(self, sphere):
        coeffR = sphere.radius * math.sqrt(3.0) / 3.0
        p1 = self.invertOnPoint(sphere.center + np.array([coeffR, coeffR, coeffR]))
        p2 = self.invertOnPoint(sphere.center + np.array([-coeffR, -coeffR, -coeffR]))
        p3 = self.invertOnPoint(sphere.center + np.array([coeffR, -coeffR, -coeffR]))
        p4 = self.invertOnPoint(sphere.center + np.array([coeffR, coeffR, -coeffR]))
        return self.__makeSphereFromPoints(p1, p2, p3, p4)


centers_out = []
radii_out = []
for base_c, base_r, centers, radii in zip(base_center, base_radius, centers_in, radii_in):
    base_sphere = Sphere(base_c[0], base_r[0])
    surroundings = []
    spheres = []
    tag = []
    for i, (c, r) in enumerate(zip(centers, radii)):
        sphere = Sphere(c, r)
        surroundings.append(sphere)

        newSphere = sphere.invertOnSphere(base_sphere)
        spheres.append(newSphere)
        tag.append(i)

    num = [0, len(spheres)]
    for level in range(0, max_level):
        for i in range(num[level], len(spheres)):
            for target_index in range(0, len(surroundings)):
                if tag[i] == target_index:
                    continue
                newSphere = surroundings[target_index].invertOnSphere(spheres[i])
                tag.append(target_index)
                spheres.append(newSphere)
                
        start_index = len(spheres)
        num.append(start_index)
        
    centers_out.append([s.center.tolist() for s in spheres])
    radii_out.append([s.radius for s in spheres])