#***************************************************************************
#*                                                                         *
#*   Copyright (c) 2009, 2010                                              *
#*   Yorik van Havre <yorik@uncreated.net>, Ken Cline <cline@frii.com>     *  
#*                                                                         *
#*   This program is free software; you can redistribute it and/or modify  *
#*   it under the terms of the GNU Lesser General Public License (LGPL)    *
#*   as published by the Free Software Foundation; either version 2 of     *
#*   the License, or (at your option) any later version.                   *
#*   for detail see the LICENCE text file.                                 *
#*                                                                         *
#*   This program is distributed in the hope that it will be useful,       *
#*   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
#*   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
#*   GNU Library General Public License for more details.                  *
#*                                                                         *
#*   You should have received a copy of the GNU Library General Public     *
#*   License along with this program; if not, write to the Free Software   *
#*   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  *
#*   USA                                                                   *
#*                                                                         *
#***************************************************************************

__title__="FreeCAD Draft Workbench - Geometry library"
__author__ = "Yorik van Havre, Jacques-Antoine Gaudin, Ken Cline"
__url__ = ["http://free-cad.sourceforge.net"]

"this file contains generic geometry functions for manipulating Part shapes"

import FreeCAD, Part, DraftVecUtils, math, cmath
from FreeCAD import Vector

NORM = Vector(0,0,1) # provisory normal direction for all geometry ops.

params = FreeCAD.ParamGet("User parameter:BaseApp/Preferences/Mod/Draft")
precision = params.GetInt("precision")

# Generic functions *********************************************************


def vec(edge):
	"vec(edge) or vec(line) -- returns a vector from an edge or a Part.line"
	# if edge is not straight, you'll get strange results!
	if isinstance(edge,Part.Shape):
		return edge.Vertexes[-1].Point.sub(edge.Vertexes[0].Point)
	elif isinstance(edge,Part.Line):
		return edge.EndPoint.sub(edge.StartPoint)
	else:
		return None

def edg(p1,p2):
	"edg(Vector,Vector) -- returns an edge from 2 vectors"
	if isinstance(p1,FreeCAD.Vector) and isinstance(p2,FreeCAD.Vector):
		if DraftVecUtils.equals(p1,p2): return None
		else: return Part.Line(p1,p2).toShape()

def getVerts(shape):
        "getVerts(shape) -- returns a list containing vectors of each vertex of the shape"
        p = []
        for v in shape.Vertexes:
                p.append(v.Point)
        return p

def v1(edge):
	"v1(edge) -- returns the first point of an edge"
	return edge.Vertexes[0].Point

def isNull(something):
        '''returns true if the given shape is null or the given placement is 0 or
        if the given vector is (0,0,0)'''
        if isinstance(something,Part.Shape):
                return something.isNull()
        elif isinstance(something,FreeCAD.Vector):
                if something == Vector(0,0,0):
                        return True
                else:
                        return False
        elif isinstance(something,FreeCAD.Placement):
                if (something.Base == Vector(0,0,0)) and (something.Rotation.Q == (0,0,0,1)):
                        return True
                else:
                        return False

def isPtOnEdge(pt,edge) :
	'''isPtOnEdge(Vector,edge) -- Tests if a point is on an edge'''
	if isinstance(edge.Curve,Part.Line) :
		orig = edge.Vertexes[0].Point
		if DraftVecUtils.isNull(pt.sub(orig).cross(vec(edge))) :
			return pt.sub(orig).Length <= vec(edge).Length and pt.sub(orig).dot(vec(edge)) >= 0
		else : 
			return False
		
	elif isinstance(edge.Curve,Part.Circle) :
		center = edge.Curve.Center
		axis   = edge.Curve.Axis ; axis.normalize()
		radius = edge.Curve.Radius
		if  round(pt.sub(center).dot(axis),precision) == 0 \
		and round(pt.sub(center).Length - radius,precision) == 0 :
			if len(edge.Vertexes) == 1 :
				return True # edge is a complete circle
			else :
				begin = edge.Vertexes[0].Point
				end   = edge.Vertexes[-1].Point
				if DraftVecUtils.isNull(pt.sub(begin)) or DraftVecUtils.isNull(pt.sub(end)) :
					return True
				else :
					# newArc = Part.Arc(begin,pt,end)
					# return DraftVecUtils.isNull(newArc.Center.sub(center)) \
					#    and DraftVecUtils.isNull(newArc.Axis-axis) \
					#    and round(newArc.Radius-radius,precision) == 0
					angle1 = DraftVecUtils.angle(begin.sub(center))
					angle2 = DraftVecUtils.angle(end.sub(center))
					anglept = DraftVecUtils.angle(pt.sub(center))
					if (angle1 < anglept) and (anglept < angle2):
						return True
	return False

def hasCurves(shape):
        "checks if the given shape has curves"
        for e in shape.Edges:
                if not isinstance(e.Curve,Part.Line):
                        return True
        return False

def isAligned(edge,axis="x"):
        "checks if the given edge or line is aligned to the given axis (x, y or z)"
        if axis == "x":
                if isinstance(edge,Part.Edge):
                        if len(edge.Vertexes) == 2:
                                if edge.Vertexes[0].X == edge.Vertexes[-1].X:
                                        return True
                elif isinstance(edge,Part.Line):
                        if edge.StartPoint.x == edge.EndPoint.x:
                                        return True
        elif axis == "y":
                if isinstance(edge,Part.Edge):
                        if len(edge.Vertexes) == 2:
                                if edge.Vertexes[0].Y == edge.Vertexes[-1].Y:
                                        return True
                elif isinstance(edge,Part.Line):
                        if edge.StartPoint.y == edge.EndPoint.y:
                                        return True
        elif axis == "z":
                if isinstance(edge,Part.Edge):
                        if len(edge.Vertexes) == 2:
                                if edge.Vertexes[0].Z == edge.Vertexes[-1].Z:
                                        return True
                elif isinstance(edge,Part.Line):
                        if edge.StartPoint.z == edge.EndPoint.z:
                                        return True
        return False

def hasOnlyWires(shape):
        "returns True if all the edges are inside a wire"
        ne = 0
        for w in shape.Wires:
                ne += len(w.Edges)
        if ne == len(shape.Edges):
                return True
        return False
                        
# edge functions *****************************************************************

def findEdge(anEdge,aList):
        '''findEdge(anEdge,aList): returns True if anEdge is found in aList of edges'''
        for e in range(len(aList)):
                if str(anEdge.Curve) == str(aList[e].Curve):
                        if DraftVecUtils.equals(anEdge.Vertexes[0].Point,aList[e].Vertexes[0].Point):
                                if DraftVecUtils.equals(anEdge.Vertexes[-1].Point,aList[e].Vertexes[-1].Point):
                                        return(e)
        return None
	

def findIntersection(edge1,edge2,infinite1=False,infinite2=False,ex1=False,ex2=False) :
	'''findIntersection(edge1,edge2,infinite1=False,infinite2=False):
	returns a list containing the intersection point(s) of 2 edges.
        You can also feed 4 points instead of edge1 and edge2'''

        pt1 = None
	
        if isinstance(edge1,FreeCAD.Vector) and isinstance(edge2,FreeCAD.Vector):
                # we got points directly
                pt1 = edge1
                pt2 = edge2
                pt3 = infinite1
                pt4 = infinite2
                infinite1 = ex1
                infinite2 = ex2

        elif isinstance(edge1.Curve,Part.Line) and isinstance(edge2.Curve,Part.Line) :	
		# we have 2 edges	
		pt1, pt2, pt3, pt4 = [edge1.Vertexes[0].Point,
                                      edge1.Vertexes[1].Point,
                                      edge2.Vertexes[0].Point,
                                      edge2.Vertexes[1].Point]

        if pt1:
                # first check if we don't already have coincident endpoints
                if (pt1 in [pt3,pt4]):
                        return [pt1]
                elif (pt2 in [pt3,pt4]):
                        return [pt2]
                
                #we have 2 straight lines		  
		if DraftVecUtils.isNull(pt2.sub(pt1).cross(pt3.sub(pt1)).cross(pt2.sub(pt4).cross(pt3.sub(pt4)))):
			vec1 = pt2.sub(pt1) ; vec2 = pt4.sub(pt3)
                        if DraftVecUtils.isNull(vec1) or DraftVecUtils.isNull(vec2):
                                return []
			vec1.normalize()  ; vec2.normalize()
			cross = vec1.cross(vec2)
			if not DraftVecUtils.isNull(cross) :
				k = ((pt3.z-pt1.z)*(vec2.x-vec2.y)+(pt3.y-pt1.y)*(vec2.z-vec2.x)+ \
					 (pt3.x-pt1.x)*(vec2.y-vec2.z))/(cross.x+cross.y+cross.z)
				vec1.scale(k,k,k)
				int = pt1.add(vec1)
				
				if  infinite1 == False and not isPtOnEdge(int,edge1) :
					return []
					
				if  infinite2 == False and not isPtOnEdge(int,edge2) :
					return []
					
				return [int]
			else :
				return [] # Lines have same direction
		else :
			return [] # Lines aren't on same plane
			
	elif isinstance(edge1.Curve,Part.Circle) and isinstance(edge2.Curve,Part.Line) \
	  or isinstance(edge1.Curve,Part.Line)   and isinstance(edge2.Curve,Part.Circle) :
	  	
	  	# deals with an arc or circle and a line
	  	
		edges = [edge1,edge2]
		for edge in edges :
			if isinstance(edge.Curve,Part.Line) :
				line = edge
			else :
				arc  = edge
				
		dirVec = vec(line) ; dirVec.normalize()
		pt1    = line.Vertexes[0].Point
		pt2    = line.Vertexes[1].Point
                pt3    = arc.Vertexes[0].Point
                pt4    = arc.Vertexes[-1].Point
		center = arc.Curve.Center

                # first check for coincident endpoints
                if (pt1 in [pt3,pt4]):
                        return [pt1]
                elif (pt2 in [pt3,pt4]):
                        return [pt2]
                
		if DraftVecUtils.isNull(pt1.sub(center).cross(pt2.sub(center)).cross(arc.Curve.Axis)) :
			# Line and Arc are on same plane
			
			dOnLine = center.sub(pt1).dot(dirVec)
			onLine  = Vector(dirVec) ; onLine.scale(dOnLine,dOnLine,dOnLine)
			toLine  = pt1.sub(center).add(onLine)
			
			if toLine.Length < arc.Curve.Radius :
				dOnLine = (arc.Curve.Radius**2 - toLine.Length**2)**(0.5)
				onLine  = Vector(dirVec) ; onLine.scale(dOnLine,dOnLine,dOnLine)
				int = [center.add(toLine).add(onLine)]
				onLine  = Vector(dirVec) ; onLine.scale(-dOnLine,-dOnLine,-dOnLine)
				int += [center.add(toLine).add(onLine)]
			elif round(toLine.Length-arc.Curve.Radius,precision) == 0 :
				int = [center.add(toLine)]
			else :
				return []
				
		else : # Line isn't on Arc's plane
			if dirVec.dot(arc.Curve.Axis) != 0 :
				toPlane  = Vector(arc.Curve.Axis) ; toPlane.normalize()
				d = pt1.dot(toPlane)
				dToPlane = center.sub(pt1).dot(toPlane)
				toPlane = Vector(pt1)
				toPlane.scale(dToPlane/d,dToPlane/d,dToPlane/d)
				ptOnPlane = toPlane.add(pt1)
				if round(ptOnPlane.sub(center).Length - arc.Curve.Radius,precision) == 0 :
					int = [ptOnPlane]
				else :
					return []
			else :
				return []
				
		if infinite1 == False :
			for i in range(len(int)-1,-1,-1) :
				if not isPtOnEdge(int[i],edge1) :
					del int[i]
		if infinite2 == False :
			for i in range(len(int)-1,-1,-1) :
				if not isPtOnEdge(int[i],edge2) :
					del int[i]
		return int
		
	elif isinstance(edge1.Curve,Part.Circle) and isinstance(edge2.Curve,Part.Circle) :
		
		# deals with 2 arcs or circles

		cent1, cent2 = edge1.Curve.Center, edge2.Curve.Center
		rad1 , rad2  = edge1.Curve.Radius, edge2.Curve.Radius
		axis1, axis2 = edge1.Curve.Axis  , edge2.Curve.Axis
		c2c          = cent2.sub(cent1)
		
		if DraftVecUtils.isNull(axis1.cross(axis2)) :
			if round(c2c.dot(axis1),precision) == 0 :
				# circles are on same plane
				dc2c = c2c.Length ;
				if not DraftVecUtils.isNull(c2c): c2c.normalize()
				if round(rad1+rad2-dc2c,precision) < 0 \
				or round(rad1-dc2c-rad2,precision) > 0 or round(rad2-dc2c-rad1,precision) > 0 :
					return []
				else :
					norm = c2c.cross(axis1)
					if not DraftVecUtils.isNull(norm): norm.normalize()
					if DraftVecUtils.isNull(norm): x = 0
					else: x = (dc2c**2 + rad1**2 - rad2**2)/(2*dc2c)
					y = abs(rad1**2 - x**2)**(0.5)
					c2c.scale(x,x,x)
					if round(y,precision) != 0 :
						norm.scale(y,y,y)
						int =  [cent1.add(c2c).add(norm)]
						int += [cent1.add(c2c).sub(norm)]
					else :
						int = [cent1.add(c2c)]
			else :
				return [] # circles are on parallel planes
		else :
			# circles aren't on same plane
			axis1.normalize() ; axis2.normalize()
			U = axis1.cross(axis2)
			V = axis1.cross(U)
			dToPlane = c2c.dot(axis2)
			d        = V.add(cent1).dot(axis2)
			V.scale(dToPlane/d,dToPlane/d,dToPlane/d)
			PtOn2Planes = V.add(cent1)
			planeIntersectionVector = U.add(PtOn2Planes)
			intTemp = findIntersection(planeIntersectionVector,edge1,True,True)
			int = []
			for pt in intTemp :
				if round(pt.sub(cent2).Length-rad2,precision) == 0 :
					int += [pt]
					
		if infinite1 == False :
			for i in range(len(int)-1,-1,-1) :
				if not isPtOnEdge(int[i],edge1) :
					del int[i]
		if infinite2 == False :
			for i in range(len(int)-1,-1,-1) :
				if not isPtOnEdge(int[i],edge2) :
					del int[i]
		
		return int
	else :
		print "fcgeo: Unsupported curve type: (" + str(edge1.Curve) + ", " + str(edge2.Curve) + ")"

def geom(edge,plac=FreeCAD.Placement()):
    "returns a Line, ArcOfCircle or Circle geom from the given edge, according to the given placement"
    if isinstance(edge.Curve,Part.Line):
            return edge.Curve
    elif isinstance(edge.Curve,Part.Circle):
        if len(edge.Vertexes) == 1:
            return Part.Circle(edge.Curve.Center,edge.Curve.Axis,edge.Curve.Radius)
        else:
            # reorienting the arc along the correct normal
            normal = plac.Rotation.multVec(FreeCAD.Vector(0,0,1))
            v1 = edge.Vertexes[0].Point
            v2 = edge.Vertexes[-1].Point
            c = edge.Curve.Center
            cu = Part.Circle(edge.Curve.Center,normal,edge.Curve.Radius)
            ref = plac.Rotation.multVec(Vector(1,0,0))
            a1 = math.pi + DraftVecUtils.angle(v1.sub(c),ref,normal)
            a2 = DraftVecUtils.angle(v2.sub(c),ref,normal)

            # direction check
            if a1 > a2:
                a1,a2 = a2,a1
                
            #print "creating sketch arc from ",cu, ", p1=",v1, " (",math.degrees(a1), "d) p2=",v2," (", math.degrees(a2),"d)"
            p= Part.ArcOfCircle(cu,a1,a2)
            return p
    else:
        return edge.Curve

def mirror (point, edge):
	"finds mirror point relative to an edge"
	normPoint = point.add(findDistance(point, edge, False))
	if normPoint:
		normPoint_point = Vector.sub(point, normPoint)
		normPoint_refl = DraftVecUtils.neg(normPoint_point)
		refl = Vector.add(normPoint, normPoint_refl)
		return refl
	else:
		return None

def findClosest(basepoint,pointslist):
	'''
	findClosest(vector,list)
	in a list of 3d points, finds the closest point to the base point.
	an index from the list is returned.
	'''
	if not pointslist: return None
	smallest = 100000
	for n in range(len(pointslist)):
		new = basepoint.sub(pointslist[n]).Length
		if new < smallest:
			smallest = new
			npoint = n
	return npoint

def concatenate(shape):
	"concatenate(shape) -- turns several faces into one"
	edges = getBoundary(shape)
	edges = sortEdges(edges)
	try:
		wire=Part.Wire(edges)
		face=Part.Face(wire)
	except:
		print "fcgeo: Couldn't join faces into one"
		return(shape)
	else:
		if not wire.isClosed():	return(wire)
		else: return(face)

def getBoundary(shape):
	"getBoundary(shape) -- this function returns the boundary edges of a group of faces"
        # make a lookup-table where we get the number of occurrences
	# to each edge in the fused face
        if isinstance(shape,list):
                shape = Part.makeCompound(shape)
	lut={}
	for f in shape.Faces:
		for e in f.Edges:
			hc= e.hashCode()
			if lut.has_key(hc): lut[hc]=lut[hc]+1
			else: lut[hc]=1
	# filter out the edges shared by more than one sub-face
	bound=[]
	for e in shape.Edges:
		if lut[e.hashCode()] == 1: bound.append(e)
	return bound

def isLine(bsp):
        "returns True if the given BSpline curve is a straight line"
        step = bsp.LastParameter/10
        b = bsp.tangent(0)
        for i in range(10):
                if bsp.tangent(i*step) != b:
                        return False
        return True

def sortEdges(lEdges, aVertex=None):
    "an alternative, more accurate version of Part.__sortEdges__"

    #There is no reason to limit this to lines only because every non-closed edge always
    #has exactly two vertices (wmayer)
    #for e in lEdges:
    #        if not isinstance(e.Curve,Part.Line):
    #                print "Warning: sortedges cannot treat wired containing curves yet."
    #                return lEdges

    def lookfor(aVertex, inEdges):
        ''' Look for (aVertex, inEdges) returns count, the position of the instance
        the position in the instance and the instance of the Edge'''
        count = 0
        linstances = [] #lists the instances of aVertex
        for i in range(len(inEdges)) :
            for j in range(2) :
                if aVertex.Point == inEdges[i].Vertexes[j-1].Point:
                        instance = inEdges[i]
                        count += 1
                        linstances += [i,j-1,instance]
        return [count]+linstances

    if (len(lEdges) < 2):
        if aVertex == None:
            return lEdges
        else:
            result = lookfor(aVertex,lEdges)
            if result[0] != 0:
                if aVertex.Point == result[3].Vertexes[0].Point:
                    return lEdges
                else:
                    if isinstance(result[3].Curve,Part.Line):
                        return [Part.Line(aVertex.Point,result[3].Vertexes[0].Point).toShape()]
                    elif isinstance(result[3].Curve,Part.Circle):
                        mp = findMidpoint(result[3])
                        return [Part.Arc(aVertex.Point,mp,result[3].Vertexes[0].Point).toShape()]
                    elif isinstance(result[3].Curve,Part.BSplineCurve):
                        if isLine(result[3].Curve):
                            return [Part.Line(aVertex.Point,result[3].Vertexes[0].Point).toShape()]
                        else:
                            return lEdges
                    else:
                        return lEdges

    olEdges = [] # ol stands for ordered list 
    if aVertex == None:
        for i in range(len(lEdges)*2) :
            if len(lEdges[i/2].Vertexes) > 1:
                result = lookfor(lEdges[i/2].Vertexes[i%2],lEdges)
                if result[0] == 1 :  # Have we found an end ?
                    olEdges = sortEdges(lEdges, result[3].Vertexes[result[2]])
                    return olEdges
        # if the wire is closed there is no end so choose 1st Vertex
        # print "closed wire, starting from ",lEdges[0].Vertexes[0].Point
        return sortEdges(lEdges, lEdges[0].Vertexes[0]) 
    else :
        #print "looking ",aVertex.Point
        result = lookfor(aVertex,lEdges)
        if result[0] != 0 :
            del lEdges[result[1]]
            next = sortEdges(lEdges, result[3].Vertexes[-((-result[2])^1)])
            #print "result ",result[3].Vertexes[0].Point,"    ",result[3].Vertexes[1].Point, " compared to ",aVertex.Point
            if aVertex.Point == result[3].Vertexes[0].Point:
                #print "keeping"
                olEdges += [result[3]] + next
            else:
                #print "inverting", result[3].Curve
                if isinstance(result[3].Curve,Part.Line):
                    newedge = Part.Line(aVertex.Point,result[3].Vertexes[0].Point).toShape()
                    olEdges += [newedge] + next
                elif isinstance(result[3].Curve,Part.Circle):
                    mp = findMidpoint(result[3])
                    newedge = Part.Arc(aVertex.Point,mp,result[3].Vertexes[0].Point).toShape()
                    olEdges += [newedge] + next
                elif isinstance(result[3].Curve,Part.BSplineCurve):
                    if isLine(result[3].Curve):
                        newedge = Part.Line(aVertex.Point,result[3].Vertexes[0].Point).toShape()
                        olEdges += [newedge] + next
                    else:
                        olEdges += [result[3]] + next
                else:
                    olEdges += [result[3]] + next                                        
            return olEdges
        else :
            return []


def findWires(edgeslist):
        '''finds connected wires in the given list of edges'''

        def touches(e1,e2):
                if len(e1.Vertexes) < 2:
                        return False
                if len(e2.Vertexes) < 2:
                        return False
                if DraftVecUtils.equals(e1.Vertexes[0].Point,e2.Vertexes[0].Point):
                        return True
                if DraftVecUtils.equals(e1.Vertexes[0].Point,e2.Vertexes[-1].Point):
                        return True
                if DraftVecUtils.equals(e1.Vertexes[-1].Point,e2.Vertexes[0].Point):
                        return True
                if DraftVecUtils.equals(e1.Vertexes[-1].Point,e2.Vertexes[-1].Point):
                        return True
                return False
        
        edges = edgeslist[:]
        wires = []
        lost = []
        while edges:
                e = edges[0]
                if not wires:
                        # create first group
                        edges.remove(e)
                        wires.append([e])
                else:
                        found = False
                        for w in wires:
                                if not found:
                                        for we in w:
                                                if touches(e,we):
                                                        edges.remove(e)
                                                        w.append(e)
                                                        found = True
                                                        break
                        if not found:
                                if e in lost:
                                        # we already tried this edge, and still nothing
                                        edges.remove(e)
                                        wires.append([e])
                                        lost = []
                                else:
                                        # put to the end of the list
                                        edges.remove(e)
                                        edges.append(e)
                                        lost.append(e)
        nwires = []
        for w in wires:
                nwires.append(Part.Wire(w))
        return nwires
                
def superWire(edgeslist,closed=False):
        '''superWire(edges,[closed]): forces a wire between edges that don't necessarily
        have coincident endpoints. If closed=True, wire will always be closed'''
        def median(v1,v2):
                vd = v2.sub(v1)
                vd.scale(.5,.5,.5)
                return v1.add(vd)
        edges = sortEdges(edgeslist)
        print edges
        newedges = []
        for i in range(len(edges)):
                curr = edges[i]
                if i == 0:
                        if closed:
                                prev = edges[-1]
                        else:
                                prev = None
                else:
                        prev = edges[i-1]
                if i == (len(edges)-1):
                        if closed:
                                next = edges[0]
                        else:
                                next = None
                else:
                        next = edges[i+1]
                print i,prev,curr,next
                if prev:
                        if curr.Vertexes[0].Point == prev.Vertexes[-1].Point:
                                p1 = curr.Vertexes[0].Point
                        else:
                                p1 = median(curr.Vertexes[0].Point,prev.Vertexes[-1].Point)
                else:
                        p1 = curr.Vertexes[0].Point
                if next:
                        if curr.Vertexes[-1].Point == next.Vertexes[0].Point:
                                p2 = next.Vertexes[0].Point
                        else:
                                p2 = median(curr.Vertexes[-1].Point,next.Vertexes[0].Point)
                else:
                        p2 = curr.Vertexes[-1].Point
                if isinstance(curr.Curve,Part.Line):
                        print "line",p1,p2
                        newedges.append(Part.Line(p1,p2).toShape())
                elif isinstance(curr.Curve,Part.Circle):
                        p3 = findMidpoint(curr)
                        print "arc",p1,p3,p2
                        newedges.append(Part.Arc(p1,p3,p2).toShape())
                else:
                        print "Cannot superWire edges that are not lines or arcs"
                        return None
        print newedges
        return Part.Wire(newedges)

def findMidpoint(edge):
	"calculates the midpoint of an edge"
	first = edge.Vertexes[0].Point
	last = edge.Vertexes[-1].Point
	if isinstance(edge.Curve,Part.Circle):
                center = edge.Curve.Center
		radius = edge.Curve.Radius
                if len(edge.Vertexes) == 1:
                        # Circle
                        dv = first.sub(center)
                        dv = DraftVecUtils.neg(dv)
                        return center.add(dv)
		axis = edge.Curve.Axis
		chord = last.sub(first)
		perp = chord.cross(axis)
		perp.normalize()
		ray = first.sub(center)
		apothem = ray.dot(perp)
		sagitta = radius - apothem
		startpoint = Vector.add(first, DraftVecUtils.scale(chord,0.5))
		endpoint = DraftVecUtils.scaleTo(perp,sagitta)
		return Vector.add(startpoint,endpoint)

	elif isinstance(edge.Curve,Part.Line):
		halfedge = DraftVecUtils.scale(last.sub(first),.5)
		return Vector.add(first,halfedge)

	else:
		return None

def complexity(obj):
	'''
	tests given object for shape complexity:
	1: line
	2: arc
	3: circle
	4: open wire with no arc
	5: closed wire
	6: wire with arcs
	7: faces
	8: faces with arcs
	'''
	shape = obj.Shape
	if shape.Faces:
		for e in shape.Edges:
			if (isinstance(e.Curve,Part.Circle)): return 8
		return 7
	if shape.Wires:
		for e in shape.Edges:
			if (isinstance(e.Curve,Part.Circle)): return 6
		for w in shape.Wires:
			if w.isClosed(): return 5
		return 4
	if (isinstance(shape.Edges[0].Curve,Part.Circle)):
		if len(shape.Vertexes) == 1:
			return 3
		return 2
	return 1

def findPerpendicular(point,edgeslist,force=None):
	'''
	findPerpendicular(vector,wire,[force]):
	finds the shortest perpendicular distance between a point and an edgeslist.
	If force is specified, only the edge[force] will be considered, and it will be
	considered infinite.
	The function will return a list	[vector_from_point_to_closest_edge,edge_index]
	or None if no perpendicular vector could be found.
	'''
        if not isinstance(edgeslist,list):
                try:
                        edgeslist = edgeslist.Edges
                except:
                        return None
	if (force == None):
		valid = None
		for edge in edgeslist:
			dist = findDistance(point,edge,strict=True)
			if dist:
				if not valid: valid = [dist,edgeslist.index(edge)]
				else:
					if (dist.Length < valid[0].Length):
						valid = [dist,edgeslist.index(edge)]
		return valid
	else:
		edge = edgeslist[force]
		dist = findDistance(point,edge)
		if dist: return [dist,force]
		else: return None
        return None

def offset(edge,vector):
	'''
	offset(edge,vector)
	returns a copy of the edge at a certain (vector) distance
	if the edge is an arc, the vector will be added at its first point
	and a complete circle will be returned
	'''
	if (not isinstance(edge,Part.Shape)) or (not isinstance(vector,FreeCAD.Vector)):
		return None
	if isinstance(edge.Curve,Part.Line):
		v1 = Vector.add(edge.Vertexes[0].Point, vector)
		v2 = Vector.add(edge.Vertexes[-1].Point, vector)
		return Part.Line(v1,v2).toShape()
	else:
		rad = edge.Vertexes[0].Point.sub(edge.Curve.Center)
		newrad = Vector.add(rad,vector).Length
		return Part.Circle(edge.Curve.Center,NORM,newrad).toShape()

def isReallyClosed(wire):
        "checks if a wire is really closed"
        if len(wire.Edges) == len(wire.Vertexes): return True
        v1 = wire.Vertexes[0].Point
        v2 = wire.Vertexes[-1].Point
        if DraftVecUtils.equals(v1,v2): return True
        return False

def getNormal(shape):
        "finds the normal of a shape, if possible"
        n = Vector(0,0,1)
        if (shape.ShapeType == "Face") and hasattr(shape,"normalAt"):
                n = shape.normalAt(0.5,0.5)
        elif shape.ShapeType == "Edge":
                if isinstance(shape.Curve,Part.Circle):
                        n = shape.Curve.Axis
        else:
                for e in shape.Edges:
                        if isinstance(e.Curve,Part.Circle):
                                n = e.Curve.Axis
                                break
                        e1 = vec(shape.Edges[0])
                        for i in range(1,len(shape.Edges)):
                                e2 = vec(shape.Edges[i])
                                if 0.1 < abs(e1.getAngle(e2)) < 1.56:
                                        n = e1.cross(e2).normalize()
                                        break
        if FreeCAD.GuiUp:
    		import Draft
	        vdir = Draft.get3DView().getViewDirection()
	        if n.getAngle(vdir) < 0.78: n = DraftVecUtils.neg(n)
        return n

def getRotation(v1,v2=FreeCAD.Vector(0,0,1)):
    '''Get the rotation Quaternion between 2 vectors'''
    if (v1.dot(v2) > 0.999999) or (v1.dot(v2) < -0.999999):
        # vectors are opposite
        return None
    axis = v1.cross(v2)
    axis.normalize()
    angle = math.degrees(math.sqrt((v1.Length ^ 2) * (v2.Length ^ 2)) + v1.dot(v2))
    return FreeCAD.Rotation(axis,angle)

def calculatePlacement(shape):
    '''calculatePlacement(shape): if the given shape is planar, this function
    returns a placement located at the center of gravity of the shape, and oriented
    towards the shape's normal. Otherwise, it returns a null placement.'''
    if not isPlanar(shape):
        return FreeCAD.Placement()
    pos = shape.BoundBox.Center
    norm = getNormal(shape)
    pla = FreeCAD.Placement()
    pla.Base = pos
    r =  getRotation(norm)
    if r:
        pla.Rotation = r
    return pla

def offsetWire(wire,dvec,bind=False,occ=False):
        '''
        offsetWire(wire,vector,[bind]): offsets the given wire along the
        given vector. The vector will be applied at the first vertex of
        the wire. If bind is True (and the shape is open), the original
        wire and the offsetted one are bound by 2 edges, forming a face.
        '''
        edges = sortEdges(wire.Edges)
        norm = getNormal(wire)
        closed = isReallyClosed(wire)
        nedges = []
        if occ:
                l=abs(dvec.Length)
                if not l: return None
                if wire.Wires:
                		wire = wire.Wires[0]
                else:
                        wire = Part.Wire(edges)
                try:
                        off = wire.makeOffset(l)
                except:
                        return None
                else:
                        return off
        for i in range(len(edges)):
                curredge = edges[i]
                delta = dvec
                if i != 0:
                        angle = DraftVecUtils.angle(vec(edges[0]),vec(curredge),norm)
                        delta = DraftVecUtils.rotate(delta,angle,norm)
                nedge = offset(curredge,delta)
                nedges.append(nedge)
        nedges = connect(nedges,closed)
        if bind and not closed:
                e1 = Part.Line(edges[0].Vertexes[0].Point,nedges[0].Vertexes[0].Point).toShape()
                e2 = Part.Line(edges[-1].Vertexes[-1].Point,nedges[-1].Vertexes[-1].Point).toShape()
                alledges = edges.extend(nedges)
                alledges = alledges.extend([e1,e2])
                w = Part.Wire(alledges)
                return w
        else:
                return nedges

def connect(edges,closed=False):
        '''connects the edges in the given list by their intersections'''
        nedges = []
        for i in range(len(edges)):
                curr = edges[i]
                # print "fcgeo.connect edge ",i," : ",curr.Vertexes[0].Point,curr.Vertexes[-1].Point
                if i > 0:
                        prev = edges[i-1]
                else:
                        if closed:
                                prev = edges[-1]
                        else:
                                prev = None
                if i < (len(edges)-1):
                        next = edges[i+1]
                else:
                        if closed: next = edges[0]
                        else:
                                next = None
                if prev:
                        # print "debug: fcgeo.connect prev : ",prev.Vertexes[0].Point,prev.Vertexes[-1].Point
                        i = findIntersection(curr,prev,True,True)
                        if i:
                                v1 = i[0]
                        else:
                                v1 = curr.Vertexes[0].Point
                else:
                        v1 = curr.Vertexes[0].Point
                if next:
                        # print "debug: fcgeo.connect next : ",next.Vertexes[0].Point,next.Vertexes[-1].Point
                        i = findIntersection(curr,next,True,True)
                        if i:
                                v2 = i[0]
                        else:
                                v2 = curr.Vertexes[-1].Point 
                else:
                        v2 = curr.Vertexes[-1].Point
                if isinstance(curr.Curve,Part.Line):
                        if v1 != v2:
                                nedges.append(Part.Line(v1,v2).toShape())
                elif isinstance(curr.Curve,Part.Circle):
                        if v1 != v2:
                                nedges.append(Part.Arc(v1,findMidPoint(curr),v2))
        try:
                return Part.Wire(nedges)
        except:
                return None

def findDistance(point,edge,strict=False):
	'''
	findDistance(vector,edge,[strict]) - Returns a vector from the point to its
	closest point on the edge. If strict is True, the vector will be returned
	only if its endpoint lies on the edge.
	'''
	if isinstance(point, FreeCAD.Vector):
		if isinstance(edge.Curve, Part.Line):
			segment = vec(edge)
			chord = edge.Vertexes[0].Point.sub(point)
			norm = segment.cross(chord)
			perp = segment.cross(norm)
			dist = DraftVecUtils.project(chord,perp)
			if not dist: return None
			newpoint = point.add(dist)
			if (dist.Length == 0):
				return None
			if strict:
				s1 = newpoint.sub(edge.Vertexes[0].Point)
				s2 = newpoint.sub(edge.Vertexes[-1].Point)
				if (s1.Length <= segment.Length) and (s2.Length <= segment.Length):
					return dist
				else:
					return None
			else: return dist
		elif isinstance(edge.Curve, Part.Circle):
			ve1 = edge.Vertexes[0].Point
			if (len(edge.Vertexes) > 1):
				ve2 = edge.Vertexes[-1].Point
			else:
				ve2 = None
			center = edge.Curve.Center
			segment = center.sub(point)
			ratio = (segment.Length - edge.Curve.Radius) / segment.Length
			dist = DraftVecUtils.scale(segment,ratio)
			newpoint = Vector.add(point, dist)
			if (dist.Length == 0):
				return None
			if strict and ve2:
				ang1 = DraftVecUtils.angle(ve1.sub(center))
				ang2 = DraftVecUtils.angle(ve2.sub(center))
				angpt = DraftVecUtils.angle(newpoint.sub(center))
				if ((angpt <= ang2 and angpt >= ang1) or (angpt <= ang1 and angpt >= ang2)):
					return dist
				else:
					return None
			else:
				return dist
                elif isinstance(edge.Curve,Part.BSplineCurve):
                        try:
                                pr = edge.Curve.parameter(point)
                                np = edge.Curve.value(pr)
                                dist = np.sub(point)
                        except:
                                print "fcgeo: Unable to get curve parameter for point ",point
                                return None
                        else:
                                return dist
		else:
			print "fcgeo: Couldn't project point"
			return None
	else:
		print "fcgeo: Couldn't project point"
		return None


def angleBisection(edge1, edge2):
	"angleBisection(edge,edge) - Returns an edge that bisects the angle between the 2 edges."
	if isinstance(edge1.Curve, Part.Line) and isinstance(edge2.Curve, Part.Line):
		p1 = edge1.Vertexes[0].Point
		p2 = edge1.Vertexes[-1].Point
		p3 = edge2.Vertexes[0].Point
		p4 = edge2.Vertexes[-1].Point
		int = findIntersection(edge1, edge2, True, True)
		if int:
			line1Dir = p2.sub(p1)
			angleDiff = DraftVecUtils.angle(line1Dir, p4.sub(p3))
			ang = angleDiff * 0.5
			origin = int[0]
			line1Dir.normalize()
			dir = DraftVecUtils.rotate(line1Dir, ang)
			return Part.Line(origin,origin.add(dir)).toShape()
		else:
			diff = p3.sub(p1)
			origin = p1.add(DraftVecUtils.scale(diff, 0.5))
			dir = p2.sub(p1); dir.normalize()
			return Part.Line(origin,origin.add(dir)).toShape()
	else:
		return None

def findClosestCircle(point,circles):
	"findClosestCircle(Vector, list of circles) -- returns the circle with closest center"
	dist = 1000000
	closest = None
	for c in circles:
		if c.Center.sub(point).Length < dist:
			dist = c.Center.sub(point).Length
			closest = c
	return closest

def isCoplanar(faces):
        "checks if all faces in the given list are coplanar"
        if len(faces) < 2:
                return True
        base =faces[0].normalAt(0,0)
        for i in range(1,len(faces)):
            for v in faces[i].Vertexes:
                chord = v.Point.sub(faces[0].Vertexes[0].Point)
                dist = DraftVecUtils.project(chord,base)
                if round(dist.Length,DraftVecUtils.precision()) > 0:
                    return False
        return True

def isPlanar(shape):
	"checks if the given shape is planar"
	if len(shape.Vertexes) <= 3:
		return True
	pts = [v.Point for v in shape.Vertexes[0:3]]
	bt = Part.Face(Part.makePolygon(pts+[pts[0]]))
	n = bt.normalAt(0,0)
	for p in shape.Vertexes[3:]:
		pv = p.Point.sub(pts[0])
		rv = DraftVecUtils.project(pv,n)
		if not DraftVecUtils.isNull(rv):
			return False
	return True

def findWiresOld(edges):
        '''finds connected edges in the list, and returns a list of lists containing edges
        that can be connected'''
        def verts(shape):
                return [shape.Vertexes[0].Point,shape.Vertexes[-1].Point]
        def group(shapes):
                shapesIn = shapes[:]
                shapesOut = [shapesIn.pop()]
                changed = False
                for s in shapesIn:
                        if len(s.Vertexes) < 2:
                                continue
                        else:
                                clean = True
                                for v in verts(s):
                                        for i in range(len(shapesOut)):
                                                if clean and (v in verts(shapesOut[i])):
                                                        shapesOut[i] = Part.Wire(shapesOut[i].Edges+s.Edges)
                                                        changed = True
                                                        clean = False
                                if clean:
                                        shapesOut.append(s)
                return(changed,shapesOut)
        working = True
        edgeSet = edges
        while working:
                result = group(edgeSet)
                working = result[0]
                edgeSet = result[1]
        return result[1]

def getTangent(edge,frompoint=None):
        '''
        returns the tangent to an edge. If from point is given, it is used to
        calculate the tangent (only useful for an arc of course).
        '''
        if isinstance(edge.Curve,Part.Line):
                return vec(edge)
        elif isinstance(edge.Curve,Part.BSplineCurve):
                if not frompoint:
                        return None
                cp = edge.Curve.parameter(frompoint)
                return edge.Curve.tangent(cp)[0]
        elif isinstance(edge.Curve,Part.Circle):
                if not frompoint:
                        v1 = edge.Vertexes[0].Point.sub(edge.Curve.Center)
                else:
                        v1 = frompoint.sub(edge.Curve.Center)
                return v1.cross(edge.Curve.Axis)
        return None

def bind(w1,w2):
        '''bind(wire1,wire2): binds 2 wires by their endpoints and
        returns a face'''
        if w1.isClosed() and w2.isClosed():
                d1 = w1.BoundBox.DiagonalLength
                d2 = w2.BoundBox.DiagonalLength
                if d1 > d2:
                        #w2.reverse()
                        return Part.Face([w1,w2])
                else:
                        #w1.reverse()
                        return Part.Face([w2,w1])
        else:
                w3 = Part.Line(w1.Vertexes[0].Point,w2.Vertexes[0].Point).toShape()
                w4 = Part.Line(w1.Vertexes[-1].Point,w2.Vertexes[-1].Point).toShape()
                return Part.Face(Part.Wire(w1.Edges+[w3]+w2.Edges+[w4]))

def cleanFaces(shape):
        "removes inner edges from coplanar faces"
        faceset = shape.Faces
        def find(hc):
                "finds a face with the given hashcode"
                for f in faceset:
                        if f.hashCode() == hc:
                                return f

        def findNeighbour(hface,hfacelist):
                "finds the first neighbour of a face in a list, and returns its index"
                eset = []
                for e in find(hface).Edges:
                        eset.append(e.hashCode())
                for i in range(len(hfacelist)):
                        for ee in find(hfacelist[i]).Edges:
                                if ee.hashCode() in eset:
                                        return i
                return None
        
        # build lookup table
        lut = {}
        for face in faceset:
                for edge in face.Edges:
                        if edge.hashCode() in lut:
                                lut[edge.hashCode()].append(face.hashCode())
                        else:
                                lut[edge.hashCode()] = [face.hashCode()]
        # print "lut:",lut
        # take edges shared by 2 faces
        sharedhedges = []
        for k,v in lut.iteritems():
                if len(v) == 2:
                        sharedhedges.append(k)
        # print len(sharedhedges)," shared edges:",sharedhedges
        # find those with same normals
        targethedges = []
        for hedge in sharedhedges:
                faces = lut[hedge]
                n1 = find(faces[0]).normalAt(0.5,0.5)
                n2 = find(faces[1]).normalAt(0.5,0.5)
                if n1 == n2:
                        targethedges.append(hedge)
        # print len(targethedges)," target edges:",targethedges
        # get target faces
        hfaces = []
        for hedge in targethedges:
                for f in lut[hedge]:
                        if not f in hfaces:
                                hfaces.append(f)

        # print len(hfaces)," target faces:",hfaces
        # sort islands
        islands = [[hfaces.pop(0)]]
        currentisle = 0
        currentface = 0
        found = True
        while hfaces:
                if not found:
                        if len(islands[currentisle]) > (currentface + 1):
                                currentface += 1
                                found = True
                        else:
                                islands.append([hfaces.pop(0)])
                                currentisle += 1
                                currentface = 0
                                found = True
                else:
                        f = findNeighbour(islands[currentisle][currentface],hfaces)
                        if f != None:
                                islands[currentisle].append(hfaces.pop(f))
                        else:
                                found = False
        # print len(islands)," islands:",islands
        # make new faces from islands
        newfaces = []
        treated = []
        for isle in islands:
                treated.extend(isle)
                fset = []
                for i in isle: fset.append(find(i))
                bounds = getBoundary(fset)
                shp = Part.Wire(sortEdges(bounds))
                shp = Part.Face(shp)
                if shp.normalAt(0.5,0.5) != find(isle[0]).normalAt(0.5,0.5):
                        shp.reverse()
                newfaces.append(shp)
        # print "new faces:",newfaces
        # add remaining faces
        for f in faceset:
                if not f.hashCode() in treated:
                        newfaces.append(f)
        # print "final faces"
        # finishing
        fshape = Part.makeShell(newfaces)
        if shape.isClosed():
                fshape = Part.makeSolid(fshape)
        return fshape


def isCubic(shape):
    '''isCubic(shape): verifies if a shape is cubic, that is, has
    8 vertices, 6 faces, and all angles are 90 degrees.'''
    # first we try fast methods
    if len(shape.Vertexes) != 8:
        return False
    if len(shape.Faces) != 6:
        return False
    if len(shape.Edges) != 12:
        return False
    for e in shape.Edges:
        if not isinstance(e.Curve,Part.Line):
            return False
    # if ok until now, let's do more advanced testing
    for f in shape.Faces:
        if len(f.Edges) != 4: return False
        for i in range(4):
            e1 = vec(f.Edges[i])
            if i < 3:
                e2 = vec(f.Edges[i+1])
            else: e2 = vec(f.Edges[0])
            rpi = [0.0,round(math.pi/2,precision)]
            if not round(e1.getAngle(e2),precision) in rpi:
                return False
    return True

def getCubicDimensions(shape):
    '''getCubicDimensions(shape): returns a list containing the placement,
    the length, the width and the height of a cubic shape. If not cubic, nothing
    is returned. The placement point is the lowest corner of the shape.'''
    if not isCubic(shape): return None
    # determine lowest face, which will be our base
    z = [10,1000000000000]
    for i in range(len(shape.Faces)):
        if shape.Faces[i].CenterOfMass.z < z[1]:
            z = [i,shape.Faces[i].CenterOfMass.z]
    if z[0] > 5: return None
    base = shape.Faces[z[0]]
    basepoint = base.Edges[0].Vertexes[0].Point
    plpoint = base.CenterOfMass
    basenorm = base.normalAt(0.5,0.5)
    # getting length and width
    vx = vec(base.Edges[0])
    vy = vec(base.Edges[1])
    # getting rotations
    rotZ = DraftVecUtils.angle(vx)
    rotY = DraftVecUtils.angle(vx,FreeCAD.Vector(vx.x,vx.y,0))
    rotX = DraftVecUtils.angle(vy,FreeCAD.Vector(vy.x,vy.y,0))
    # getting height
    vz = None
    rpi = round(math.pi/2,precision)
    for i in range(1,6):
        for e in shape.Faces[i].Edges:
            if basepoint in [e.Vertexes[0].Point,e.Vertexes[1].Point]:
                vtemp = vec(e)
                # print vtemp
                if round(vtemp.getAngle(vx),precision) == rpi:
                    if round(vtemp.getAngle(vy),precision) == rpi:
                        vz = vtemp
    if not vz: return None
    mat = FreeCAD.Matrix()
    mat.move(plpoint)
    mat.rotateX(rotX)
    mat.rotateY(rotY)
    mat.rotateZ(rotZ)
    return [FreeCAD.Placement(mat),round(vx.Length,precision),round(vy.Length,precision),round(vz.Length,precision)]

def removeInterVertices(wire):
        '''removeInterVertices(wire) - remove unneeded vertices (those that
        are in the middle of a straight line) from a wire, returns a new wire.'''
        edges = sortEdges(wire.Edges)
        nverts = []
        def getvec(v1,v2):
                if not abs(round(v1.getAngle(v2),precision) in [0,round(math.pi,precision)]):
                        nverts.append(edges[i].Vertexes[-1].Point)
        for i in range(len(edges)-1):
                vA = vec(edges[i])
                vB = vec(edges[i+1])
                getvec(vA,vB)
        vA = vec(edges[-1])
        vB = vec(edges[0])
        getvec(vA,vB)
        if nverts:
                if wire.isClosed():
                        nverts.append(nverts[0])
                w = Part.makePolygon(nverts)
                return w
        else:
                return wire

def arcFromSpline(edge):
        """arcFromSpline(edge): turns the given edge into an arc, by taking
        its first point, midpoint and endpoint. Works best with bspline
        segments such as those from imported svg files. Use this only
        if you are sure your edge is really an arc..."""
        if isinstance(edge.Curve,Part.Line):
                print "This edge is straight, cannot build an arc on it"
                return None
        if len(edge.Vertexes) > 1:
                # 2-point arc
                p1 = edge.Vertexes[0].Point
                p2 = edge.Vertexes[-1].Point
                ml = edge.Length/2
                p3 = edge.valueAt(ml)
                try:
                        return Part.Arc(p1,p3,p2).toShape()
                except:
                        print "Couldn't make an arc out of this edge"
                        return None
        else:
                # circle
                p1 = edge.Vertexes[0].Point
                ml = edge.Length/2
                p2 = edge.valueAt(ml)
                ray = p2.sub(p1)
                ray.scale(.5,.5,.5)
                center = p1.add(ray)
                radius = ray.Length
                try:
                        return Part.makeCircle(radius,center)
                except:
                        print "couldn't make a circle out of this edge"

# Fillet code graciously donated by Jacques-Antoine Gaudin

def fillet(lEdges,r):
	''' Take a list of two Edges & a float as argument,
		Returns a list of sorted edges describing a round corner'''

        def getCurveType(edge,existingCurveType = None):
                '''Builds or completes a dictionnary containing edges with keys "Arc" and "Line"'''	
                if not existingCurveType :
                        existingCurveType = { 'Line' : [], 'Arc' : [] }
                if issubclass(type(edge.Curve),Part.Line) :
                        existingCurveType['Line'] += [edge]
                elif issubclass(type(edge.Curve),Part.Circle) :
                        existingCurveType['Arc']  += [edge]
                else :
                        raise Exception("Edge's curve must be either Line or Arc")
                return existingCurveType
	
	rndEdges = lEdges[0:2]
	rndEdges = sortEdges(rndEdges)

	if len(rndEdges) < 2 :
		return rndEdges
	
	if r <= 0 :
		print "fcgeo.fillet : Error : radius is negative."
		return rndEdges
		
	curveType = getCurveType(rndEdges[0])
	curveType = getCurveType(rndEdges[1],curveType)
		
	lVertexes = rndEdges[0].Vertexes + [rndEdges[1].Vertexes[-1]]
		
	if len(curveType['Line']) == 2:
		
		# Deals with 2-line-edges lists -------------------------------------- 
		
		U1 = lVertexes[0].Point.sub(lVertexes[1].Point) ; U1.normalize()
		U2 = lVertexes[2].Point.sub(lVertexes[1].Point) ; U2.normalize()
		alpha = U1.getAngle(U2)
		
		if round(alpha,precision) == 0 or round(alpha - math.pi,precision) == 0: # Edges have same direction
			print "fcgeo.fillet : Warning : edges have same direction. Did nothing"
			return rndEdges
			
		dToCenter = r / math.sin(alpha/2.) 
		dToTangent = (dToCenter**2-r**2)**(0.5)
		dirVect = Vector(U1) ; dirVect.scale(dToTangent,dToTangent,dToTangent)
		arcPt1 = lVertexes[1].Point.add(dirVect)
		
		dirVect = U2.add(U1) ; dirVect.normalize()
		dirVect.scale(dToCenter-r,dToCenter-r,dToCenter-r)
		arcPt2 = lVertexes[1].Point.add(dirVect)
		
		dirVect = Vector(U2) ; dirVect.scale(dToTangent,dToTangent,dToTangent)
		arcPt3 = lVertexes[1].Point.add(dirVect)
			
		if (dToTangent>lEdges[0].Length) or (dToTangent>lEdges[1].Length) :
			print "fcgeo.fillet : Error : radius value ", r," is too high"
			return rndEdges
		
		rndEdges[1]   =  Part.Edge(Part.Arc(arcPt1,arcPt2,arcPt3))
		rndEdges[0]   =  Part.Edge(Part.Line(lVertexes[0].Point,arcPt1))
		rndEdges     += [Part.Edge(Part.Line(arcPt3,lVertexes[2].Point))]
		
		return rndEdges
		
	elif len(curveType['Arc']) == 1 :
		
		# Deals with lists containing an arc and a line ----------------------------------
		
		if lEdges[0] in curveType['Arc'] :
			lineEnd = lVertexes[2] ; arcEnd = lVertexes[0] ; arcFirst = True
		else :
			lineEnd = lVertexes[0] ; arcEnd = lVertexes[2] ; arcFirst = False
		arcCenter = curveType['Arc'][0].Curve.Center
		arcRadius = curveType['Arc'][0].Curve.Radius
		arcAxis   = curveType['Arc'][0].Curve.Axis
		arcLength = curveType['Arc'][0].Length
		
		U1 = lineEnd.Point.sub(lVertexes[1].Point) ; U1.normalize()
		toCenter = arcCenter.sub(lVertexes[1].Point)
		if arcFirst : # make sure the tangent points towards the arc
			T = arcAxis.cross(toCenter)
		else :
			T = toCenter.cross(arcAxis)
			
		projCenter = toCenter.dot(U1)
		if round(abs(projCenter),precision) > 0 :
			normToLine = U1.cross(T).cross(U1)
		else :
			normToLine = Vector(toCenter) 
		normToLine.normalize()
		
		dCenterToLine = toCenter.dot(normToLine) - r
		
		if  round(projCenter,precision) > 0 :
			newRadius = arcRadius - r
		elif round(projCenter,precision) < 0 or (round(projCenter,precision) == 0 and U1.dot(T) > 0):
			newRadius = arcRadius + r
		else :
			print "fcgeo.fillet : Warning : edges are already tangent. Did nothing"
			return rndEdges
		
		toNewCent = newRadius**2-dCenterToLine**2
		if toNewCent > 0 :
			toNewCent = abs(abs(projCenter) - toNewCent**(0.5))
		else :
			print "fcgeo.fillet : Error : radius value ", r," is too high"
			return rndEdges
			
		U1.scale(toNewCent,toNewCent,toNewCent)
		normToLine.scale(r,r,r)	
		newCent = lVertexes[1].Point.add(U1).add(normToLine)
		
		arcPt1= lVertexes[1].Point.add(U1)
		arcPt2= lVertexes[1].Point.sub(newCent); arcPt2.normalize()
		arcPt2.scale(r,r,r) ; arcPt2 = arcPt2.add(newCent)
		if newRadius == arcRadius - r :
			arcPt3= newCent.sub(arcCenter)
		else :
			arcPt3= arcCenter.sub(newCent) 
		arcPt3.normalize()
		arcPt3.scale(r,r,r) ; arcPt3 = arcPt3.add(newCent)
		arcPt = [arcPt1,arcPt2,arcPt3]
		
		
		# Warning : In the following I used a trick for calling the right element
		# in arcPt or V : arcFirst is a boolean so - not arcFirst is -0 or -1
		# list[-1] is the last element of a list and list[0] the first
		# this way I don't have to proceed tests to know the position of the arc
		
		myTrick = not arcFirst
		
		V  = [arcPt3]
		V += [arcEnd.Point]
		
		toCenter.scale(-1,-1,-1)
		
		delLength = arcRadius * V[0].sub(arcCenter).getAngle(toCenter)
		if delLength > arcLength or toNewCent > curveType['Line'][0].Length: 
			print "fcgeo.fillet : Error : radius value ", r," is too high"
			return rndEdges
			
		arcAsEdge = arcFrom2Pts(V[-arcFirst],V[-myTrick],arcCenter,arcAxis)
		
		V = [lineEnd.Point,arcPt1]
		lineAsEdge = Part.Edge(Part.Line(V[-arcFirst],V[myTrick]))
		
		rndEdges[not arcFirst]   =  arcAsEdge
		rndEdges[arcFirst]       =  lineAsEdge
		rndEdges[1:1] = [Part.Edge(Part.Arc(arcPt[- arcFirst],arcPt[1],arcPt[- myTrick]))]
		
		return rndEdges
		
	elif len(curveType['Arc']) == 2 :
		
		# Deals with lists of 2 arc-edges --------------------------------------------
		
		arcCenter, arcRadius, arcAxis, arcLength, toCenter, T, newRadius = [], [], [], [], [], [], []
		for i in range(2) :
			arcCenter += [curveType['Arc'][i].Curve.Center]
			arcRadius += [curveType['Arc'][i].Curve.Radius]
			arcAxis   += [curveType['Arc'][i].Curve.Axis]
			arcLength += [curveType['Arc'][i].Length]
			toCenter  += [arcCenter[i].sub(lVertexes[1].Point)]
		T += [arcAxis[0].cross(toCenter[0])]
		T += [toCenter[1].cross(arcAxis[1])]
		CentToCent = toCenter[1].sub(toCenter[0])
		dCentToCent = CentToCent.Length
		
		sameDirection = (arcAxis[0].dot(arcAxis[1]) > 0)
		TcrossT = T[0].cross(T[1])
		if sameDirection :
			if   round(TcrossT.dot(arcAxis[0]),precision) > 0 :
				newRadius += [arcRadius[0]+r]
				newRadius += [arcRadius[1]+r]
			elif round(TcrossT.dot(arcAxis[0]),precision) < 0 :
				newRadius += [arcRadius[0]-r]
				newRadius += [arcRadius[1]-r]
			elif T[0].dot(T[1]) > 0 :
				newRadius += [arcRadius[0]+r]
				newRadius += [arcRadius[1]+r]
			else :
				print "fcgeo.fillet : Warning : edges are already tangent. Did nothing"
				return rndEdges
		elif not sameDirection : 
			if   round(TcrossT.dot(arcAxis[0]),precision) > 0 :
				newRadius += [arcRadius[0]+r]
				newRadius += [arcRadius[1]-r]
			elif round(TcrossT.dot(arcAxis[0]),precision) < 0 :
				newRadius += [arcRadius[0]-r]
				newRadius += [arcRadius[1]+r]
			elif T[0].dot(T[1]) > 0 :
				if arcRadius[0] > arcRadius[1] :
					newRadius += [arcRadius[0]-r]
					newRadius += [arcRadius[1]+r]
				elif arcRadius[1] > arcRadius[0] :
					newRadius += [arcRadius[0]+r]
					newRadius += [arcRadius[1]-r]
				else :
					print "fcgeo.fillet : Warning : arcs are coincident. Did nothing"
					return rndEdges
			else :
				print "fcgeo.fillet : Warning : edges are already tangent. Did nothing"
				return rndEdges
				
		if newRadius[0]+newRadius[1] < dCentToCent or \
		   newRadius[0]-newRadius[1] > dCentToCent or \
		   newRadius[1]-newRadius[0] > dCentToCent :
		   	print "fcgeo.fillet : Error : radius value ", r," is too high"
			return rndEdges
				
		x = (dCentToCent**2+newRadius[0]**2-newRadius[1]**2)/(2*dCentToCent)
		y = (newRadius[0]**2-x**2)**(0.5)
		
		CentToCent.normalize() ; toCenter[0].normalize() ; toCenter[1].normalize()
		if abs(toCenter[0].dot(toCenter[1])) != 1 :
			normVect = CentToCent.cross(CentToCent.cross(toCenter[0]))
		else :
			normVect = T[0]
		normVect.normalize()
		CentToCent.scale(x,x,x) ; normVect.scale(y,y,y)
		newCent = arcCenter[0].add(CentToCent.add(normVect))
		CentToNewCent = [newCent.sub(arcCenter[0]),newCent.sub(arcCenter[1])]
		for i in range(2) :
			CentToNewCent[i].normalize()
			if newRadius[i] == arcRadius[i]+r :
				CentToNewCent[i].scale(-r,-r,-r)
			else :
				CentToNewCent[i].scale(r,r,r)
		toThirdPt = lVertexes[1].Point.sub(newCent) ; toThirdPt.normalize()
		toThirdPt.scale(r,r,r)
		arcPt1 = newCent.add(CentToNewCent[0])
		arcPt2 = newCent.add(toThirdPt)
		arcPt3 = newCent.add(CentToNewCent[1])
		arcPt = [arcPt1,arcPt2,arcPt3]
		
		arcAsEdge = []
		for i in range(2) :
			toCenter[i].scale(-1,-1,-1)
			delLength = arcRadius[i] * arcPt[-i].sub(arcCenter[i]).getAngle(toCenter[i])
			if delLength > arcLength[i] :
				print "fcgeo.fillet : Error : radius value ", r," is too high"
				return rndEdges
			V = [arcPt[-i],lVertexes[-i].Point]
			arcAsEdge += [arcFrom2Pts(V[i-1],V[-i],arcCenter[i],arcAxis[i])]
		
		rndEdges[0]   =  arcAsEdge[0]
		rndEdges[1]   =  arcAsEdge[1]
		rndEdges[1:1] = [Part.Edge(Part.Arc(arcPt[0],arcPt[1],arcPt[2]))]
		
		return rndEdges
	       		
def filletWire(aWire,r,makeClosed=True):
	''' Fillets each angle of a wire with r as radius value'''
	
	edges = aWire.Edges
	edges = sortEdges(edges)
	filEdges = [edges[0]]
	for i in range(len(edges)-1):
		result = fillet([filEdges[-1],edges[i+1]],r)
		if len(result)>2:
			filEdges[-1:] = result[0:3]
		else :
			filEdges[-1:] = result[0:2]
	if isReallyClosed(aWire) and makeClosed :
		result = fillet([filEdges[-1],filEdges[0]],r)
		if len(result)>2:
			filEdges[-1:] = result[0:2]
			filEdges[0]   = result[2]
	return Part.Wire(filEdges)

   
# circle functions *********************************************************

def getBoundaryAngles(angle,alist):
        '''returns the 2 closest angles from the list that
        encompass the given angle'''
        negs = True
        while negs:
                negs = False
                for i in range(len(alist)):
                        if alist[i] < 0:
                                alist[i] = 2*math.pi + alist[i]
                                negs = True
                if angle < 0:
                        angle = 2*math.pi + angle
                        negs = True
        lower = None
        for a in alist:
                if a < angle:
                        if lower == None:
                                lower = a
                        else:
                                if a > lower:
                                        lower = a
        if lower == None:
                lower = 0
                for a in alist:
                        if a > lower:
                                lower = a
        higher = None
        for a in alist:
                if a > angle:
                        if higher == None:
                                higher = a
                        else:
                                if a < higher:
                                        higher = a
        if higher == None:
                higher = 2*math.pi
                for a in alist:
                        if a < higher:
                                higher = a
        return (lower,higher)
                        

def circleFrom2tan1pt(tan1, tan2, point):
	"circleFrom2tan1pt(edge, edge, Vector)"
	if isinstance(tan1.Curve, Part.Line) and isinstance(tan2.Curve, Part.Line) and isinstance(point, FreeCAD.Vector):
		return circlefrom2Lines1Point(tan1, tan2, point)
	elif isinstance(tan1.Curve, Part.Circle) and isinstance(tan2.Curve, Part.Line) and isinstance(point, FreeCAD.Vector):
		return circlefromCircleLinePoint(tan1, tan2, point)
	elif isinstance(tan2.Curve, Part.Circle) and isinstance(tan1.Curve, Part.Line) and isinstance(point, FreeCAD.Vector):
		return circlefromCircleLinePoint(tan2, tan1, point)
	elif isinstance(tan2.Curve, Part.Circle) and isinstance(tan1.Curve, Part.Circle) and isinstance(point, FreeCAD.Vector):
		return circlefrom2Circles1Point(tan2, tan1, point)

def circleFrom2tan1rad(tan1, tan2, rad):
	"circleFrom2tan1rad(edge, edge, float)"
	if isinstance(tan1.Curve, Part.Line) and isinstance(tan2.Curve, Part.Line):
		return circleFrom2LinesRadius(tan1, tan2, rad)
	elif isinstance(tan1.Curve, Part.Circle) and isinstance(tan2.Curve, Part.Line):
		return circleFromCircleLineRadius(tan1, tan2, rad)
	elif isinstance(tan1.Curve, Part.Line) and isinstance(tan2.Curve, Part.Circle):
		return circleFromCircleLineRadius(tan2, tan1, rad)
	elif isinstance(tan1.Curve, Part.Circle) and isinstance(tan2.Curve, Part.Circle):
		return circleFrom2CirclesRadius(tan1, tan2, rad)

def circleFrom1tan2pt(tan1, p1, p2):
	if isinstance(tan1.Curve, Part.Line) and isinstance(p1, FreeCAD.Vector) and isinstance(p2, FreeCAD.Vector):
		return circlefrom1Line2Points(tan1, p1, p2)
	if isinstance(tan1.Curve, Part.Line) and isinstance(p1, FreeCAD.Vector) and isinstance(p2, FreeCAD.Vector):
		return circlefrom1Circle2Points(tan1, p1, p2)

def circleFrom1tan1pt1rad(tan1, p1, rad):
	if isinstance(tan1.Curve, Part.Line) and isinstance(p1, FreeCAD.Vector):
		return circleFromPointLineRadius(p1, tan1, rad)
	if isinstance(tan1.Curve, Part.Circle) and isinstance(p1, FreeCAD.Vector):
		return circleFromPointCircleRadius(p1, tan1, rad)

def circleFrom3tan(tan1, tan2, tan3):
	tan1IsLine = isinstance(tan1.Curve, Part.Line)
	tan2IsLine = isinstance(tan2.Curve, Part.Line)
	tan3IsLine = isinstance(tan3.Curve, Part.Line)
	tan1IsCircle = isinstance(tan1.Curve, Part.Circle)
	tan2IsCircle = isinstance(tan2.Curve, Part.Circle)
	tan3IsCircle = isinstance(tan3.Curve, Part.Circle)
	if tan1IsLine and tan2IsLine and tan3IsLine:
		return circleFrom3LineTangents(tan1, tan2, tan3)
	elif tan1IsCircle and tan2IsCircle and tan3IsCircle:
		return circleFrom3CircleTangents(tan1, tan2, tan3)
	elif (tan1IsCircle and tan2IsLine and tan3IsLine):
		return circleFrom1Circle2Lines(tan1, tan2, tan3)
	elif (tan1IsLine and tan2IsCircle and tan3IsLine):
		return circleFrom1Circle2Lines(tan2, tan1, tan3)
	elif (tan1IsLine and tan2IsLine and tan3IsCircle):
		return circleFrom1Circle2Lines(tan3, tan1, tan2)
	elif (tan1IsLine and tan2IsCircle and tan3IsCircle):
		return circleFrom2Circle1Lines(tan2, tan3, tan1)
	elif (tan1IsCircle and tan2IsLine and tan3IsCircle):
		return circleFrom2Circle1Lines(tan1, tan3, tan2)
	elif (tan1IsCircle and tan2IsCircle and tan3IsLine):
		return circleFrom2Circle1Lines(tan1, tan2, tan3)

def circlefrom2Lines1Point(edge1, edge2, point):
	"circlefrom2Lines1Point(edge, edge, Vector)"
	bis = angleBisection(edge1, edge2)
	if not bis: return None
	mirrPoint = mirror(point, bis)
	return circlefrom1Line2Points(edge1, point, mirrPoint)

def circlefrom1Line2Points(edge, p1, p2):
	"circlefrom1Line2Points(edge, Vector, Vector)"
	p1_p2 = edg(p1, p2)
	s = findIntersection(edge, p1_p2, True, True)
	if not s: return None
	s = s[0]
	v1 = p1.sub(s)
	v2 = p2.sub(s)
	projectedDist = math.sqrt(abs(v1.dot(v2)))
	edgeDir = vec(edge); edgeDir.normalize()
	projectedCen1 = Vector.add(s, DraftVecUtils.scale(edgeDir, projectedDist))
	projectedCen2 = Vector.add(s, DraftVecUtils.scale(edgeDir, -projectedDist))
	perpEdgeDir = edgeDir.cross(Vector(0,0,1))
	perpCen1 = Vector.add(projectedCen1, perpEdgeDir)
	perpCen2 = Vector.add(projectedCen2, perpEdgeDir)
	mid = findMidpoint(p1_p2)
	x = DraftVecUtils.crossproduct(vec(p1_p2)); x.normalize()
	perp_mid = Vector.add(mid, x)
	cen1 = findIntersection(edg(projectedCen1, perpCen1), edg(mid, perp_mid), True, True)
	cen2 = findIntersection(edg(projectedCen2, perpCen2), edg(mid, perp_mid), True, True)
	circles = []
	if cen1:
		radius = DraftVecUtils.dist(projectedCen1, cen1[0])
		circles.append(Part.Circle(cen1[0], NORM, radius))
	if cen2:
		radius = DraftVecUtils.dist(projectedCen2, cen2[0])
		circles.append(Part.Circle(cen2[0], NORM, radius))

	if circles: return circles
	else: return None

def circleFrom2LinesRadius (edge1, edge2, radius):
	"circleFrom2LinesRadius(edge,edge,radius)"
	int = findIntersection(edge1, edge2, True, True)
	if not int: return None
	int = int[0]
	bis12 = angleBisection(edge1,edge2)
	bis21 = Part.Line(bis12.Vertexes[0].Point,DraftVecUtils.rotate(vec(bis12), math.pi/2.0))
	ang12 = abs(DraftVecUtils.angle(vec(edge1),vec(edge2)))
	ang21 = math.pi - ang12
	dist12 = radius / math.sin(ang12 * 0.5)
	dist21 = radius / math.sin(ang21 * 0.5)
	circles = []
	cen = Vector.add(int, DraftVecUtils.scale(vec(bis12), dist12))
	circles.append(Part.Circle(cen, NORM, radius))
	cen = Vector.add(int, DraftVecUtils.scale(vec(bis12), -dist12))
	circles.append(Part.Circle(cen, NORM, radius))
	cen = Vector.add(int, DraftVecUtils.scale(vec(bis21), dist21))
	circles.append(Part.Circle(cen, NORM, radius))
	cen = Vector.add(int, DraftVecUtils.scale(vec(bis21), -dist21))
	circles.append(Part.Circle(cen, NORM, radius))
	return circles

def circleFrom3LineTangents (edge1, edge2, edge3):
	"circleFrom3LineTangents(edge,edge,edge)"
	def rot(ed):
		return Part.Line(v1(ed),v1(ed).add(DraftVecUtils.rotate(vec(ed),math.pi/2))).toShape()
	bis12 = angleBisection(edge1,edge2)
	bis23 = angleBisection(edge2,edge3)
	bis31 = angleBisection(edge3,edge1)
	intersections = []
	int = findIntersection(bis12, bis23, True, True)	
	if int:
		radius = findDistance(int[0],edge1).Length
		intersections.append(Part.Circle(int[0],NORM,radius))
	int = findIntersection(bis23, bis31, True, True)
	if int:
		radius = findDistance(int[0],edge1).Length
		intersections.append(Part.Circle(int[0],NORM,radius))
	int = findIntersection(bis31, bis12, True, True)
	if int:
		radius = findDistance(int[0],edge1).Length
		intersections.append(Part.Circle(int[0],NORM,radius))
	int = findIntersection(rot(bis12), rot(bis23), True, True)
	if int:
		radius = findDistance(int[0],edge1).Length
		intersections.append(Part.Circle(int[0],NORM,radius))
	int = findIntersection(rot(bis23), rot(bis31), True, True)
	if int:
		radius = findDistance(int[0],edge1).Length
		intersections.append(Part.Circle(int[0],NORM,radius))
	int = findIntersection(rot(bis31), rot(bis12), True, True)
	if int:
		radius = findDistance(int[0],edge1).Length
		intersections.append(Part.Circle(int[0],NORM,radius))
	circles = []
	for int in intersections:
		exists = False
		for cir in circles:
			if DraftVecUtils.equals(cir.Center, int.Center):
				exists = True
				break
		if not exists:
			circles.append(int)
	if circles:
		return circles
	else:
		return None

def circleFromPointLineRadius (point, edge, radius):
 	"circleFromPointLineRadius (point, edge, radius)"
	dist = findDistance(point, edge, False)
	center1 = None
	center2 = None
	if dist.Length == 0:
		segment = vec(edge)
		perpVec = DraftVecUtils.crossproduct(segment); perpVec.normalize()
		normPoint_c1 = DraftVecUtils.scale(perpVec, radius)
		normPoint_c2 = DraftVecUtils.scale(perpVec, -radius)
		center1 = point.add(normPoint_c1)
		center2 = point.add(normPoint_c2)
	elif dist.Length > 2 * radius:
		return None
	elif dist.Length == 2 * radius:
		normPoint = point.add(findDistance(point, edge, False))
		dummy = DraftVecUtils.scale(normPoint.sub(point), 0.5)
		cen = point.add(dummy)
		circ = Part.Circle(cen, NORM, radius)
		if circ:
			return [circ]
		else:
			return None
	else:
		normPoint = point.add(findDistance(point, edge, False))
		normDist = DraftVecUtils.dist(normPoint, point)
		dist = math.sqrt(radius**2 - (radius - normDist)**2)
		centerNormVec = DraftVecUtils.scaleTo(point.sub(normPoint), radius)
		edgeDir = edge.Vertexes[0].Point.sub(normPoint); edgeDir.normalize()
		center1 = centerNormVec.add(normPoint.add(DraftVecUtils.scale(edgeDir, dist)))
		center2 = centerNormVec.add(normPoint.add(DraftVecUtils.scale(edgeDir, -dist)))
	circles = []
	if center1:
		circ = Part.Circle(center1, NORM, radius)
		if circ:
			circles.append(circ)
	if center2:
		circ = Part.Circle(center2, NORM, radius)
		if circ:
			circles.append(circ)

	if len(circles):
		return circles
	else:
		return None

def circleFrom2PointsRadius(p1, p2, radius):
	"circleFrom2PointsRadiust(Vector, Vector, radius)"
	if DraftVecUtils.equals(p1, p2): return None

	p1_p2 = Part.Line(p1, p2).toShape()
	dist_p1p2 = DraftVecUtils.dist(p1, p1)
	mid = findMidpoint(p1_p2)
	if dist_p1p2 == 2*radius:
		circle = Part.Circle(mid, norm, radius)
		if circle: return [circle]
		else: return None
	dir = vec(p1_p2); dir.normalize()
	perpDir = dir.cross(Vector(0,0,1)); perpDir.normailze()
	dist = math.sqrt(radius**2 - (dist_p1p2 / 2.0)**2)
	cen1 = Vector.add(mid, DraftVecUtils.scale(perpDir, dist))
	cen2 = Vector.add(mid, DraftVecUtils.scale(perpDir, -dist))
	circles = []
	if cen1: circles.append(Part.Circle(cen1, norm, radius))
	if cen2: circles.append(Part.Circle(cen2, norm, radius))
	if circles: return circles
	else: return None



#############################33 to include







def outerSoddyCircle(circle1, circle2, circle3):
	'''
	Computes the outer soddy circle for three tightly packed circles.
	'''
	if isinstance(circle1.Curve, Part.Circle) and isinstance(circle2.Curve, Part.Circle) and isinstance(circle3.Curve, Part.Circle):
		# Original Java code Copyright (rc) 2008 Werner Randelshofer
		# Converted to python by Martin Buerbaum 2009
		# http://www.randelshofer.ch/treeviz/
		# Either Creative Commons Attribution 3.0, the MIT license, or the GNU Lesser General License LGPL.

		A = circle1.Curve.Center
		B = circle2.Curve.Center
		C = circle3.Curve.Center

		ra = circle1.Curve.Radius
		rb = circle2.Curve.Radius
		rc = circle3.Curve.Radius

		# Solution using Descartes' theorem, as described here:
		# http://en.wikipedia.org/wiki/Descartes%27_theorem
		k1 = 1 / ra
		k2 = 1 / rb
		k3 = 1 / rc
		k4 = abs(k1 + k2 + k3 - 2 * math.sqrt(k1 * k2 + k2 * k3 + k3 * k1))

		q1 = (k1 + 0j) * (A.x + A.y * 1j)
		q2 = (k2 + 0j) * (B.x + B.y * 1j)
		q3 = (k3 + 0j) * (C.x + C.y * 1j)

		temp = ((q1 * q2) + (q2 * q3) + (q3 * q1))
		q4 = q1 + q2 + q3 - ((2 + 0j) * cmath.sqrt(temp) )

		z = q4 / (k4 + 0j)

		# If the formula is not solveable, we return no circle.
		if (not z or not (1 / k4)):
			return None

		X = -z.real
		Y = -z.imag
		print "Outer Soddy circle: " + str(X) + " " + str(Y) + "\n"	# Debug

		# The Radius of the outer soddy circle can also be calculated with the following formula:
		# radiusOuter = abs(r1*r2*r3 / (r1*r2 + r1*r3 + r2*r3 - 2 * math.sqrt(r1*r2*r3 * (r1+r2+r3))))
		circ = Part.Circle(Vector(X, Y, A.z), norm, 1 / k4)
		return circ

	else:
		print "debug: outerSoddyCircle bad parameters!\n"
		# FreeCAD.Console.PrintMessage("debug: outerSoddyCircle bad parameters!\n")
		return None

def innerSoddyCircle(circle1, circle2, circle3):
	'''
	Computes the inner soddy circle for three tightly packed circles.
	'''
	if isinstance(circle1.Curve, Part.Circle) and isinstance(circle2.Curve, Part.Circle) and isinstance(circle3.Curve, Part.Circle):
		# Original Java code Copyright (rc) 2008 Werner Randelshofer
		# Converted to python by Martin Buerbaum 2009
		# http://www.randelshofer.ch/treeviz/

		A = circle1.Curve.Center
		B = circle2.Curve.Center
		C = circle3.Curve.Center

		ra = circle1.Curve.Radius
		rb = circle2.Curve.Radius
		rc = circle3.Curve.Radius

		# Solution using Descartes' theorem, as described here:
		# http://en.wikipedia.org/wiki/Descartes%27_theorem
		k1 = 1 / ra
		k2 = 1 / rb
		k3 = 1 / rc
		k4 = abs(k1 + k2 + k3 + 2 * math.sqrt(k1 * k2 + k2 * k3 + k3 * k1))

		q1 = (k1 + 0j) * (A.x + A.y * 1j)
		q2 = (k2 + 0j) * (B.x + B.y * 1j)
		q3 = (k3 + 0j) * (C.x + C.y * 1j)

		temp = ((q1 * q2) + (q2 * q3) + (q3 * q1))
		q4 = q1 + q2 + q3 + ((2 + 0j) * cmath.sqrt(temp) )

		z = q4 / (k4 + 0j)

		# If the formula is not solveable, we return no circle.
		if (not z or not (1 / k4)):
			return None

		X = z.real
		Y = z.imag
		print "Outer Soddy circle: " + str(X) + " " + str(Y) + "\n"	# Debug

		# The Radius of the inner soddy circle can also be calculated with the following formula:
		# radiusInner = abs(r1*r2*r3 / (r1*r2 + r1*r3 + r2*r3 + 2 * math.sqrt(r1*r2*r3 * (r1+r2+r3))))
		circ = Part.Circle(Vector(X, Y, A.z), norm, 1 / k4)
		return circ

	else:
		print "debug: innerSoddyCircle bad parameters!\n"
		# FreeCAD.Console.PrintMessage("debug: innerSoddyCircle bad parameters!\n")
		return None
 
def circleFrom3CircleTangents(circle1, circle2, circle3):
	'''
	http://en.wikipedia.org/wiki/Problem_of_Apollonius#Inversive_methods
	http://mathworld.wolfram.com/ApolloniusCircle.html
	http://mathworld.wolfram.com/ApolloniusProblem.html
	'''

	if isinstance(circle1.Curve, Part.Circle) and isinstance(circle2.Curve, Part.Circle) and isinstance(circle3.Curve, Part.Circle):
		int12 = findIntersection(circle1, circle2, True, True)
		int23 = findIntersection(circle2, circle3, True, True)
		int31 = findIntersection(circle3, circle1, True, True)

		if int12 and int23 and int31:
			if len(int12) == 1 and len(int23) == 1 and len(int31) == 1:
				# Only one intersection with each circle.
				# => "Soddy Circle" - 2 solutions.
				# http://en.wikipedia.org/wiki/Problem_of_Apollonius#Mutually_tangent_given_circles:_Soddy.27s_circles_and_Descartes.27_theorem
				# http://mathworld.wolfram.com/SoddyCircles.html
				# http://mathworld.wolfram.com/InnerSoddyCenter.html
				# http://mathworld.wolfram.com/OuterSoddyCenter.html

				r1 = circle1.Curve.Radius
				r2 = circle2.Curve.Radius
				r3 = circle3.Curve.Radius
				outerSoddy = outerSoddyCircle(circle1, circle2, circle3)
#				print str(outerSoddy) + "\n" # Debug

				innerSoddy = innerSoddyCircle(circle1, circle2, circle3)
#				print str(innerSoddy) + "\n" # Debug

				circles = []
				if outerSoddy:
					circles.append(outerSoddy)
				if innerSoddy:
					circles.append(innerSoddy)
				return circles

			# @todo Calc all 6 homothetic centers.
			# @todo Create 3 lines from the inner and 4 from the outer h. center.
			# @todo Calc. the 4 inversion poles of these lines for each circle.
			# @todo Calc. the radical center of the 3 circles.
			# @todo Calc. the intersection points (max. 8) of 4 lines (trough each inversion pole and the radical center) with the circle.
			#       This gives us all the tangent points.
		else:
			# Some circles are inside each other or an error has occured.
			return None

	else:
		print "debug: circleFrom3CircleTangents bad parameters!\n"
		# FreeCAD.Console.PrintMessage("debug: circleFrom3CircleTangents bad parameters!\n")
		return None


def linearFromPoints (p1, p2):
	'''
	Calculate linear equation from points.
	Calculate the slope and offset parameters of the linear equation of a line defined by two points.

	Linear equation:
	y = m * x + b
	m = dy / dx
	m ... Slope
	b ... Offset (point where the line intersects the y axis)
	dx/dy ... Delta x and y. Using both as a vector results in a non-offset direction vector.
	'''
	if isinstance(p1, Vector) and isinstance(p2, Vector):
		line = {}
		line['dx'] = (p2.x - p1.x)
		line['dy'] = (p2.y - p1.y)
		line['slope'] = line['dy'] / line['dx']
		line['offset'] = p1.y - slope * p1.x
		return line
	else:
		return None


def determinant (mat,n):
	'''
	determinant(matrix,int) - Determinat function. Returns the determinant
	of a n-matrix. It recursively expands the minors.
	'''
	matTemp = [[0.0,0.0,0.0],[0.0,0.0,0.0],[0.0,0.0,0.0]]
	if (n > 1):
		if n == 2:
			d = mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1]
		else:
			d = 0.0
			for j1 in range(n):
				# Create minor
				for i in range(1, n):
					j2 = 0
					for j in range(n):
						if j == j1:
							continue
						matTemp[i-1][j2] = mat[i][j]
						j2 += 1
				d += (-1.0)**(1.0 + j1 + 1.0) * mat[0][j1] * determinant(matTemp, n-1)
		return d
	else:
		return 0

	
def findHomotheticCenterOfCircles(circle1, circle2):
	'''
	findHomotheticCenterOfCircles(circle1, circle2)
	Calculates the homothetic center(s) of two circles.

	http://en.wikipedia.org/wiki/Homothetic_center
	http://mathworld.wolfram.com/HomotheticCenter.html
	'''

	if isinstance(circle1.Curve, Part.Circle) and isinstance(circle2.Curve, Part.Circle):
		if DraftVecUtils.equals(circle1.Curve.Center, circle2.Curve.Center):
			return None

		cen1_cen2 = Part.Line(circle1.Curve.Center, circle2.Curve.Center).toShape()
		cenDir = vec(cen1_cen2); cenDir.normalize()

		# Get the perpedicular vector.
		perpCenDir = cenDir.cross(Vector(0,0,1)); perpCenDir.normalize()

		# Get point on first circle
		p1 = Vector.add(circle1.Curve.Center, DraftVecUtils.scale(perpCenDir, circle1.Curve.Radius))

		centers = []
		# Calculate inner homothetic center
		# Get point on second circle
		p2_inner = Vector.add(circle1.Curve.Center, DraftVecUtils.scale(perpCenDir, -circle1.Curve.Radius))
		hCenterInner = DraftVecUtils.intersect(circle1.Curve.Center, circle2.Curve.Center, p1, p2_inner, True, True)
		if hCenterInner:
			centers.append(hCenterInner)

		# Calculate outer homothetic center (only exists of the circles have different radii)
		if circle1.Curve.Radius != circle2.Curve.Radius:
			# Get point on second circle
			p2_outer = Vector.add(circle1.Curve.Center, DraftVecUtils.scale(perpCenDir, circle1.Curve.Radius))
			hCenterOuter = DraftVecUtils.intersect(circle1.Curve.Center, circle2.Curve.Center, p1, p2_outer, True, True)
			if hCenterOuter:
				centers.append(hCenterOuter)

		if len(centers):
			return centers
		else:
			return None

	else:
		print "debug: findHomotheticCenterOfCircles bad parameters!\n"
		FreeCAD.Console.PrintMessage("debug: findHomotheticCenterOfCirclescleFrom3tan bad parameters!\n")
		return None


def findRadicalAxis(circle1, circle2):
	'''
	Calculates the radical axis of two circles.
	On the radical axis (also called power line) of two circles any
	tangents drawn from a point on the axis to both circles have the same length.

	http://en.wikipedia.org/wiki/Radical_axis
	http://mathworld.wolfram.com/RadicalLine.html

	@sa findRadicalCenter
	'''

	if isinstance(circle1.Curve, Part.Circle) and isinstance(circle2.Curve, Part.Circle):
		if DraftVecUtils.equals(circle1.Curve.Center, circle2.Curve.Center):
			return None
		r1 = circle1.Curve.Radius
		r2 = circle1.Curve.Radius
		cen1 = circle1.Curve.Center
		# dist .. the distance from cen1 to cen2.
		dist = DraftVecUtils.dist(cen1, circle2.Curve.Center)
		cenDir = cen1.sub(circle2.Curve.Center); cenDir.normalize()

		# Get the perpedicular vector.
		perpCenDir = cenDir.cross(Vector(0,0,1)); perpCenDir.normalize()

		# J ... The radical center.
		# K ... The point where the cadical axis crosses the line of cen1->cen2.
		# k1 ... Distance from cen1 to K.
		# k2 ... Distance from cen2 to K.
		# dist = k1 + k2

		k1 = (dist + (r1^2 - r2^2) / dist) / 2.0
		#k2 = dist - k1

		K = Vector.add(cen1, DraftVecUtils.scale(cenDir, k1))

		# K_ .. A point somewhere between K and J (actually with a distance of 1 unit from K).
		K_ = Vector,add(K, perpCenDir)

		radicalAxis = Part.Line(K, Vector.add(origin, dir))

		if radicalAxis:
			return radicalAxis
		else:
			return None
	else:
		print "debug: findRadicalAxis bad parameters!\n"
		FreeCAD.Console.PrintMessage("debug: findRadicalAxis bad parameters!\n")
		return None



def findRadicalCenter(circle1, circle2, circle3):
	'''
	findRadicalCenter(circle1, circle2, circle3):
	Calculates the radical center (also called the power center) of three circles.
	It is the intersection point of the three radical axes of the pairs of circles.

	http://en.wikipedia.org/wiki/Power_center_(geometry)
	http://mathworld.wolfram.com/RadicalCenter.html

	@sa findRadicalAxis
	'''

	if isinstance(circle1.Curve, Part.Circle) and isinstance(circle2.Curve, Part.Circle) and isinstance(circle2.Curve, Part.Circle):
		radicalAxis12 = findRadicalAxis(circle1, circle2)
		radicalAxis23 = findRadicalAxis(circle1, circle2)

		if not radicalAxis12 or not radicalAxis23:
			# No radical center could be calculated.
			return None

		int = findIntersection(radicalAxis12, radicalAxis23, True, True)

		if int:
			return int
		else:
			# No radical center could be calculated.
			return None
	else:
		print "debug: findRadicalCenter bad parameters!\n"
		FreeCAD.Console.PrintMessage("debug: findRadicalCenter bad parameters!\n")
		return None

def pointInversion(circle, point):
	'''
	pointInversion(Circle, Vector)

	Circle inversion of a point.
	Will calculate the inversed point an return it.
	If the given point is equal to the center of the circle "None" will be returned.

	See also:
	http://en.wikipedia.org/wiki/Inversive_geometry
	'''

	if isinstance(circle.Curve, Part.Circle) and isinstance(point, FreeCAD.Vector):
		cen = circle.Curve.Center
		rad = circle.Curve.Radius

		if DraftVecUtils.equals(cen, point):
			return None

		# Inverse the distance of the point
		# dist(cen -> P) = r^2 / dist(cen -> invP)

		dist = DraftVecUtils.dist(point, cen)
		invDist = rad**2 / d

		invPoint = Vector(0, 0, point.z)
		invPoint.x = cen.x + (point.x - cen.x) * invDist / dist;
		invPoint.y = cen.y + (point.y - cen.y) * invDist / dist;

		return invPoint

	else:
		print "debug: pointInversion bad parameters!\n"
		FreeCAD.Console.PrintMessage("debug: pointInversion bad parameters!\n")
		return None

def polarInversion(circle, edge):
	'''
	polarInversion(circle, edge):
	Returns the inversion pole of a line.
	edge ... The polar.
	i.e. The nearest point on the line is inversed.

	http://mathworld.wolfram.com/InversionPole.html
	'''

	if isinstance(circle.Curve, Part.Circle) and isinstance(edge.Curve, Part.Line):
		nearest = circle.Curve.Center.add(findDistance(circle.Curve.Center, edge, False))
		if nearest:
			inversionPole = pointInversion(circle, nearest)
			if inversionPole:
				return inversionPole

	else:
		print "debug: circleInversionPole bad parameters!\n"
		FreeCAD.Console.PrintMessage("debug: circleInversionPole bad parameters!\n")
		return None

def circleInversion(circle, circle2):
	'''
	pointInversion(Circle, Circle)

	Circle inversion of a circle.
	'''
	if isinstance(circle.Curve, Part.Circle) and isinstance(circle2.Curve, Part.Circle):
		cen1 = circle.Curve.Center
		rad1 = circle.Curve.Radius

		if DraftVecUtils.equals(cen1, point):
			return None

		invCen2 = Inversion(circle, circle2.Curve.Center)

		pointOnCircle2 = Vector.add(circle2.Curve.Center, Vector(circle2.Curve.Radius, 0, 0))
		invPointOnCircle2 = Inversion(circle, pointOnCircle2)

		return Part.Circle(invCen2, norm, DraftVecUtils.dist(invCen2, invPointOnCircle2))

	else:
		print "debug: circleInversion bad parameters!\n"
		FreeCAD.Console.PrintMessage("debug: circleInversion bad parameters!\n")
		return None

