__author__ = 'thoth'

#
# for
# http://blender.stackexchange.com/questions/28239/how-can-i-determine-calculate-a-circumcenter-with-3-points/28260#28260
#

import bpy
import bmesh
from mathutils import *

def replace_col(M, i, C):
    for r in range(len(M)):
        M[r][i] = C[r]

def compute_circle_center_(B, C):
    """
    Compute the center of a circle that includes points [0,0,0], B, and C
    """
    N = B.cross(C)

    m_d = Matrix( [
        B,C,N
    ])

    col = [B.dot(B)*0.5,
    C.dot(C)*0.5,
    0]

    m_x = m_d.copy()
    replace_col(m_x, 0, col)

    m_y = m_d.copy()
    replace_col(m_y, 1, col)

    m_z = m_d.copy()
    replace_col(m_z, 2, col)

    m_d_d = m_d.determinant()

    x = m_x.determinant() / m_d_d
    y = m_y.determinant() / m_d_d
    z = m_z.determinant() / m_d_d

    #print(m_d)
    #print(m_x)
    #print(m_y)
    #print(m_z)

    return Vector([x,y,z])

    latex = """\vec x \cdot \frac{\vec B}{\left | \vec B \right |} = \frac{1}{2}\left | \vec B \right |
\vec x \cdot \vec B = \frac{1}{2}\left | \vec B \right |^2
\vec x \cdot \vec C = \frac{1}{2}\left | \vec C \right |^2

\vec N = \vec B \times \vec C
\vec x \cdot \vec N = 0


x B_x + y B_y + z B_z = \frac{1}{2}(B_x^2+B_y^2+B_z^2) = b
x C_x + y C_y + z C_z = \frac{1}{2}(C_x^2+C_y^2+C_z^2) = c
x N_x + y N_y + z N_z = 0

x= \frac{\begin{vmatrix}
 b & B_y & B_z\\
 c & C_y & C_z\\
 0 & N_y & N_z
 \end{vmatrix}}{\begin{vmatrix}
 B_x & B_y & B_z\\
 C_x & C_y & C_z\\
 N_x & N_y & N_z
 \end{vmatrix}}

y= \frac{\begin{vmatrix}
 B_x & b & B_z\\
 C_x & c & C_z\\
 N_x & 0 & N_z
 \end{vmatrix}}{\begin{vmatrix}
 B_x & B_y & B_z\\
 C_x & C_y & C_z\\
 N_x & N_y & N_z
 \end{vmatrix}}

z= \frac{\begin{vmatrix}
 B_x & B_y & b\\
 C_x & C_y & c\\
 N_x & N_y & 0
 \end{vmatrix}}{\begin{vmatrix}
 B_x & B_y & B_z\\
 C_x & C_y & C_z\\
 N_x & N_y & N_z
 \end{vmatrix}}
"""

def compute_circle_center(A,B,C):
    #print ("%r\t %r\t %r"%(A,B,C))
    B_ = B-A
    C_ = C-A
    #print ("%r\t %r\t %r"%(B_,C_, N))

    return A+compute_circle_center_(B_, C_)

def space_3d():
    for area in bpy.context.screen.areas:
        for space in area.spaces:
            if space.type=='VIEW_3D':
                return space
    raise AssertionError("could not find any editor space of type VIEW_3D")

def mission():
    obj = bpy.context.active_object
    verts=[]
    if bpy.context.mode == 'EDIT_MESH':
        bm = bmesh.from_edit_mesh(obj.data)
        for v in bm.verts:
            if v.select:
                verts.append(obj.matrix_world * v.co)
    else:
        for v in obj.data.vertices:
            if v.select:
                verts.append(obj.matrix_world * v.co)
    if len(verts) != 3:
        raise AssertionError("You must select exactly 3 vertices from the mesh, not %d"%len(verts))

    space_3d().cursor_location = compute_circle_center(verts[0], verts[1], verts[2])

#
#
#

mission()