#!BPY

"""
Name: 'MegaBool'
Blender: 234
Group: 'Object'
Tooltip: 'Perform boolean operations on meshes and get nice results'
"""
__bpydoc__ ="""\
                                                                 <br>
 Experimental New Mesh Boolean Operations                        <br>
                                                                 <br>
 (C) 2004 Theodore K Schundler (tschundler)(@)(scu)(edu)         <br>
 FOR TESTING ONLY, not for redistribution                        <br>
 Lastest version available for from:                             <br>
 http://astro.scu.edu/~ted/oss/blender/boolean                   <br>
                                                                 <br>
 USAGE:                                                          <br>
   -By Menu-                                                     <br>
     * Put this file in ~/.blender/scripts/                      <br>
     * Select two objects, then choose "MegaBool" from the       <br>
        scripts menu.                                            <br>
   -By Running the script-                                       <br>
     * Load this file into blender's internal text editor        <br>
     * Select two objects, then press Alt+P in the test editor   <br>
                                                                 <br>
 OPERATIONS:                                                     <br>
   On execution of the script, there are four operations         <br>
     * Intersect - Intersection of the two meshes (logical AND)  <br>
     * Union - Union of the two meshes (logical OR)              <br>
     * Difference - First selected mesh - Last selected mesh     <br>
     * Cookie Cutter - First mesh's faces are subdivided where   <br>
        they intersect with the second mesh                      <br>
   After the script executes, all "inside" vertices are selected <br>
                                                                 <br>
 IMPORTANT NOTES:                                                <br>
   Like any boolean algorithm, normals on the source meshes must <br>
 be all facing outside. To fix normals, select a mesh, enter     <br>
 edit mode, select all verticies, and press Ctrl+N. Also, all    <br>
 intersections must create closed loops. (Basically the meshes   <br>
 should be closed.) If you find a set of meshes that do not work <br>
 send me the .blend file and an explaination of the problem.     <br>
                                                                 <br>
 TO DO:                                                          <br>
  * Set materials, UVs, etc., on new mesh                        <br>
  * Intersection of curves? (use splines for smooth shapes)      <br>

-----------------------------------------------------------------
                                                                 
  Date:        Updates:                                          <br>
  16-Aug-2004  Initial Revision - working(?) & ready to go       <br>
  17-Aug-2004  Fixed some bugs based on incorrect assumptions.   <br>
                that means even more computation. Also made more <br>
                use of the traversalCookies to speed up a some   <br>
                sections                                         <br>
  24-Aug-2004  Smart filling & other stuff...                    <br>
  26-Aug-2004  Much cleanup, nothing new / fancy                 <br>
  18-Sep-2004  Re-did filling. 2nd object can now be a set of    <br>
                loops. Other minor refinements.                  <br>
  24-Sep-2004  Fixed a wrong assumption about vertex ordering    <br>
  26-Sep-2004  Got rid of special filling stuff that made face   <br>
                selection work better because with the new edit  <br>
                modes, it isn't necissary.                       <br>
  27-Oct-2004  Made some headway with the coplanar problem, but  <br>
                there is still a colinear problem, and sometimes <br>
                an issue if a point lies in another objects plane<br>
  01-Nov-2004  Fixed edge crossing tests in joinconcentric       <br>
  10-Dec-2004  Rewrite of almost everything                      <br>
  08-Apr-2005  Rewrite of alot - completely changed how inside   <br>
                vs outside determination works                   <br>
  24-Apr-2005  Rewrite filling (handles complicated cases faster)<br>
                fixed dumb mistakes in intersecting edges, fixed <br>
                transform of vertex normals.                     <br>
                                                                 
"""

__author__ = "Theodore K. Schundler"
__version__ = "r11 24-Apr-05"
__email__ = ('Author EMail, tschundler:scu*edu')
__url__ = ('elysiun','Script Web Page, http://astro.scu.edu/~ted/oss/blender/boolean/')


import Blender
from Blender import NMesh
from Blender import Object
from Blender import Draw
from Blender import Mathutils
from Blender.Mathutils import *
import time

traversalCookie=0
messages=""
enorm=0#1
drawxpt=0#1
xpt=[]
						
###################################################################
# Utility Functions                                               #
###################################################################
def MidPoint(v1,v2):
	return (v1*0.49+v2*0.51); #a little bit of jitter ...

#Expand Bounding Box
def ExpandBounds2D(Bounds,v):
	for i in range(2):
		if (Bounds[0][i]>v[i]):
			 Bounds[0][i]=v[i]
		if (Bounds[1][i]<v[i]):
			 Bounds[1][i]=v[i]
	return Bounds

#Build a bounding box
def BuildBounds2D(v):
	Bounds=[[v[0],v[1]],[v[0],v[1]]]
	return Bounds

#Test if bounding boxes overlap
def TestBounds2D(b1,b2):
	for i in range(2):
		if (b1[0][i]-0.001>b2[1][i]) or (b2[0][i]-0.001>b1[1][i]):
			return 0
	return 1

def CrossVec2D(v):
	return Vector(-v[1],v[0])

def TestEdgeIntersectEdge2D(pA,pB,pC,pD):
	vDC=pD-pC
	vBA=pB-pA
	Base=DotVecs(vBA,CrossVec2D(vDC))
	if Base==0:
		#Parallel
		return False
	CvAC=CrossVec2D(pA-pC)
	uA=DotVecs(vDC,CvAC)/Base
	if uA>=0 and uA<=1:
		uB=DotVecs(vBA,CvAC)/Base
		if uB>=0 and uB<=1:
			return True
	return False

#Determine if an angle should be regarded as > 180 degrees
def Over180(VectorA,VectorB,EdgeVectorA,EdgeVectorB):
	#return DotVecs(VectorA+VectorB,EdgeVectorA+EdgeVectorB)<0
	return DotVecs(VectorA,EdgeVectorB)<0

def BuildBasis(normal,edge):
	e2=CrossVecs(normal,edge)
	return [edge,e2]

def Make2D(basis,vector):
	return Vector(DotVecs(basis[0],vector),DotVecs(basis[1],vector))


		
def ExpandBounds(Bounds,v):
	for i in range(3):
		if (Bounds[0][i]>v[i]):
			 Bounds[0][i]=v[i]
		if (Bounds[1][i]<v[i]):
			 Bounds[1][i]=v[i]
	return Bounds

def BuildBounds(v1,v2):
	Bounds=[[v1[0],v1[1],v1[2]],[v1[0],v1[1],v1[2]]]
	return ExpandBounds(Bounds,v2)

def TestBounds(b1,b2):
	#return 1
	for i in range(3):
		if (b1[0][i]-0.001>b2[1][i]) or (b2[0][i]-0.001>b1[1][i]):
			return 0
	return 1

def FindIntersectLineLine3D(p1,p2,p3,p4,stop_at_a=0,limits=0,advanced=0,tolerance=0.0002):
	if limits:
		#First test for shared points
		if (p1==p3 or p1==p4 or p2==p3 or p2==p4):
			return []

		#Bounds test...maybe it should be tested to see if this actually helps?
		if (TestBounds(BuildBounds(p1,p2),BuildBounds(p3,p4))==0):
			return []

	#based on http://astronomy.swin.edu.au/~pbourke/geometry/lineline3d/
	d2121=DotVecs(p2-p1,p2-p1) #note: square of the edge's length (maybe should be precomputed for all edges?)
	d4321=DotVecs(p4-p3,p2-p1)
	d4343=DotVecs(p4-p3,p4-p3) #note: square of edge's length again
	d1343=DotVecs(p1-p3,p4-p3)
	d1321=DotVecs(p1-p3,p2-p1)
	mnad = (d2121*d4343-d4321*d4321)
	if mnad==0:#mnad>-0.00000000000001 and mnad<0.00000000000001: #parallel
		#print "Parallel!!!",mnad
		return []
	else:
		mna = (d1343*d4321-d1321*d4343)/mnad
		if (limits==0 or (mna>-0.000001 and mna <1.000001)):
			mnb = (d1343+mna*d4321)/d4343
			if ((limits==0 and stop_at_a==0) or (mnb > -0.00001 and mnb <1.000001)):
				if mnad<0.001 or advanced:
					dv=p1+(p2-p1)*mna-(p3+(p4-p3)*mnb)
					dis=DotVecs(dv,dv)
					#Needs to be 0.00001 for animtest to work...
					#Maybe is I break Quads into tris, this won't be a problem?
					if dis>tolerance:
						#print "Parallel - late detection",dv,dis,(p3+(p4-p3)*mnb)
						return []
					return [mna,mnb,dis]
				return [mna,mnb]
	return []

def FindDistPointToLine3D(p1,p2,p3):
	Va=p2-p1
	Vb=p3-p1
	d=DotVecs(Va,Va)
	if d!=0:
		p4= ((DotVecs(Vb,Va)/d)*Va)+p1
	else:
		p4=p1
	Vc=p4-p3
	return DotVecs(Vc,Vc)

def TestVertInEdge3D(v,e):
	vectA=e.v[1].tco-e.v[0].tco
	dot = DotVecs(vectA,vectA)
	vectB=v.tco-e.v[0].tco
	dv=DotVecs(vectA,vectB)
	if (dv>=-0.0000001 and dv<=dot+0.0000001):
		vectA.normalize()
		vectB.normalize()
		dv2=DotVecs(vectA,vectB)
		if (dv2>0.99999999 and dv2<1.00000001):
			return 1
	return 0

def TestIntersectLineLine3D(p1,p2,p3,p4):
	return len(FindIntersectLineLine3D(p1,p2,p3,p4,0,1))==2

def TestEdgeIntersectLoop3D(p0,p1,loop):
	lastv=loop[len(loop)-1]
	for v in loop:
		if TestIntersectLineLine3D(p0.tco,p1.tco,lastv.tco,v.tco):
			return 1
		lastv=v
	return 0

def TestPointSideOfLine3D(point,v0,v1,normal):
	#Note: Maybe I should store the edge's normal on the edge?
	#This function gets called from different loops, so it may be worth it
	cp=CrossVecs(normal,v1-v0)

	s=DotVecs(point-v0,cp)
	if s < 0:
		return 0
	else:
		return 1

#test if a point is within the edges of a convex face
#Basically we make sure the point is on the "inside" side
#of all the edges of the shape.
#it may be worthwhile to optimize this since uit gets a called a lot when filling
def TestPointInConvexLoop(p,l,normal):
	lastv=l[len(l)-1]
	for v in l:
		if not TestPointSideOfLine3D(p,lastv.tco,v.tco,normal):
			return 0
		lastv=v
	return 1
	
def FindIntersectRayPlane(r_point,r_direction,p_point,p_normal,positive_only=0):
	d=DotVecs(p_normal,r_direction)
	if abs(d)<0.001:
		return 0
	n=DotVecs(p_normal,p_point-r_point)
	nd=n/d
	if nd<0 and positive_only:
		return 0
	return nd

def TestPointInPlane(point,normal,facevert):
	va=point-facevert
	#normal's length should always be 1....why isn't it?
	d=DotVecs(normal,va)/normal.length
	#must be <0.001 ....apparently much less than 0.001 is a problem too
	#>=0.001 needed for animtest frame 82
	if d<0.001 and d > -0.001:
		print "Point in Plane:",d,point,normal,facevert
		return 1
	return 0
					
def TestPointInFace(p,f):
	vcount=len(f.v)
	if vcount:
		return TestPointInConvexLoop(p,f.v,f.normal)
	else: #This is an infinitely large face, so we test between the two edges
		v0=f.edges[0].v[0]
		v1=f.edges[1].v[0]
		d0=DotVecs(v1.tco-v0.tco,p-v0.tco)
		d1=DotVecs(v0.tco-v1.tco,p-v1.tco)
		if d0>0 and d1>0:
			#print "PIF",d0,d1,v0.tco,v1.tco,p
			return 1
		else:
			#print "PNIF"
			return 0

def TestIntersectRayFace(r_point,r_direction,face,positive_only=0):
	d=FindIntersectRayPlane(r_point,r_direction,face.e[0].v[0].tco,face.normal,positive_only)
	if d==0:
		return 0
	p=r_point+d*r_direction
	return TestPointInFace(p,face)

def FindIntersectEdgeFace(edge,face):
	r_point=edge.v[0].tco
	r_direction=edge.v[1].tco-r_point
	d=FindIntersectRayPlane(r_point,r_direction,face.e[0].v[0].tco,face.normal)
	if d==0:
		return 0
	#May need to fudge these valus for fuzzy intersections
	if (not edge.infinite) and (d<0 or d>1):
		return 0
	p=r_point+d*r_direction
	if TestPointInFace(p,face):
		return p
	return 0

def SuperMakeLoop(start,vA,offset,stop,depth=0):
	if depth>999:
		print "ML(): loop too deep (maybe not a bug if the mesh is complicated)"
		return []
	global traversalCookie
	start.travCookie=traversalCookie
	if start.v[0]==vA:
		vB=start.v[1]
	else:
		vB=start.v[0]
	eNormal=start.f[offset][1]
	face=start.f[offset][2]
	side=start.f[offset][0]
	print "ML():Making loop...",vA.tco,vB.tco,start,stop
	print "ML():\t",start.f
	bestE=[]
	bestAngle=-99
	for e in vB.e:
		if e!=start:
			if e==stop:
				#all done!
				print "ML(): Loop Done!"
				start.f[offset][3]=1
				return [side,[[start,vA]]]
			if e.v[0]==vB:
				vC=e.v[1]
			else:
				vC=e.v[0]
			if vC==vB:
				print "ML(): BUG: zero-length edge"
			dA=(vB.tco-vA.tco)
			dB=(vB.tco-vC.tco)
			if (dB.length<0.00001):
				print "ML(): BUG - zero length edge, separate verts - should have caught this earlier"
				MergeVert(vB,vC)
			#is there a way to do this w/o normalization?
			dA.normalize()
			dB.normalize()
			vsum=(dA+dB)/2
			ang=DotVecs(dA,dB)
			if ang>0.9999: #this is a test for parallel edges -- see stuff in FindIntersectLineLine3D to avoid sqrt
				print "ML(): BUG? Edges folded over..."
				print "ML(): \t- ",start,start.v[0].tco,start.v[1].tco
				print "ML(): \t- ",e,e.v[0].tco,e.v[1].tco
				ang=-9999999
			print "ML(): Option? ",vB.tco,"->",vC.tco,"(",ang,")"
			#Maybe add 180 to angle if new vect :dot: current vect's normal <0, then just find smallest angle?
			if ang>bestAngle:
				for fi in range(len(e.f)):
					f=e.f[fi]
					if f[2]==face and f[3]==0:
						print "ML(): Have a candidate...",vB.tco,"->",vC.tco,"(",ang,")"
						nsum=eNormal+f[1]
						nsl=nsum.length
						nsum.normalize()
						dp=DotVecs(nsum,vsum)
						dQ=DotVecs(nsum,dA)
						dW=DotVecs(nsum,dB)
						dE=DotVecs(vsum,eNormal)
						dF=DotVecs(vsum,f[1])
						#note: the nsl>0.0001 thing is for dealing with edges where ang is 180 degrees (-1)
						print "ML():\t",ang,nsl,dQ,dW,dp,dE,dF
						#if (ang>0 or nsl>0.0001) and (abs(dQ-dW)<0.001) and ang>bestAngle:
						if (ang>0 or nsl>0.0001) and ((dE>0 and dF>0) or (dE<0 and dF<0)) and ang>bestAngle:
							#print "ML():\t\t",ang,dp,dQ,dW
							bestE=[e,fi,f[0],vC]
							bestAngle=ang
	if bestE and bestE[2]!=6:
		print "ML(): found next step",bestE[0],bestE[1],bestE[3].tco
		start.f[offset][3]=1
		r=SuperMakeLoop(bestE[0],vB,bestE[1],stop,depth+1)
		if len(r)>1:
			if r[0]!=0 and r[0]!=5:
				if side!=0 and side!=5 and side!=r[0]:
					print "ML(): BUG??? State changed midloop from ",side,"to",r[0],"at",vB.tco
				if (side==0 or side==5):
					side=r[0]
			#if start.broken==0:
			#	start.setSides(side)
			#if side!=0:
			#	start.f[offset][0]=side
			lp=r[1]
			lp.append([start,vA])
			return [side,lp]
	start.f[offset][3]=0
	print "ML(): loop failed!"
	return []
				
			
###################################################################
# Pretty Quad Filler v2                                           #
###################################################################

#Objects
class LoopElement:
	def __init__(self,co,link):
		self.co=co
		self.link=link
		self.bb=0
		self.next=0
	def copy(self):
		ne=LoopElement(self.co,self.link)
		ne.next=self.next
		return ne
	def SetNext(self,next):
		self.next=next
	def GetNext(self):
		return self.next
	def Reverse(self,first,prev):
		next=self.next
		self.next=prev
		if next==first:
			next.next=self
			return
		next.Reverse(first,self)
	def GetCo(self):
		return self.co
	def GetNextCo(self):
		return self.next.co
	def GetVector(self):
		return self.next.co-self.co
	def GetEdgeNormal(self):
		return CrossVec2D(self.GetVector())
	def TestEdgeIntersect(self,ptA,ptB,bb,stop):
		if not self.bb:
			self.bb=BuildBounds2D(self.co)
			self.bb=ExpandBounds2D(self.bb,self.next.co)
		if TestBounds2D(self.bb,bb):
			if TestEdgeIntersectEdge2D(ptA,ptB,self.co,self.next.co):
				return True
		if self==stop:
			return False
		return self.next.TestEdgeIntersect(ptA,ptB,bb,stop)
	def nFill(self,n,tolerance,stop=0):
		#print "nFill(",n,stop,")"
		if stop==self:
			return False
		if not stop:
			stop=self
		c=self
		ok=1
		#test if we can make a convex loop
		
		for i in range(n-2):
			no=c.GetEdgeNormal()
			no.normalize()
			#print "\t",c.next.co,c.next.GetVector(),no
			#this can probably be optimized to jump past the convex point
			#maybe convex points can be cached?
			if DotVecs(c.GetEdgeNormal(),c.next.GetVector())<tolerance:
				#print "Convex fail at ",c.next.co
				return self.next.nFill(n,tolerance,stop)
			prev=c
			c=c.next
		#print "\t",c.next.co,self.co-c.next.co,c.GetEdgeNormal()
		if DotVecs(c.GetEdgeNormal(),self.co-c.next.co)<tolerance:
			#print "Convex fail at 2p ",c.next.co
			return self.next.nFill(n,tolerance,stop)
		#print "\t",self.co,c.next.co-self.co,self.GetEdgeNormal()
		if DotVecs(self.GetEdgeNormal(),c.next.co-self.co)<tolerance:
			#print "Convex fail at 3p ",self.co
			return self.next.nFill(n,tolerance,stop)
		lastpt=c.next
		#print "\tconvex OK"
		
		#Does it break out of the loop (if the loop size is > n)
		if self!=lastpt.next:
			if (not Over180(c.next.GetVector(),c.GetVector()*-1.0,c.next.GetEdgeNormal(),c.GetEdgeNormal())) and DotVecs(self.co-c.next.co,c.next.GetEdgeNormal())<tolerance:
				#print "Breaks out of original poly (test 1)",self.co,c.next.co,c.next.GetEdgeNormal()
				return self.next.nFill(n,tolerance,stop)
		
		#ok, so it's convex. But does it intersect any existing edges?
		c=c.next
		if c!=self:
			prev=c
			c=c.next
			#print "cross?",prev.co,c.co
			#test through all the other points in the loop to make sure
			#our face doesn't enclose them (which would mean edges intersect somewhere)
			while c!=self:
				#print "\tcross?",c.co
				tp=self
				ok=0
				#Note: should probably use and epsilon instead of zero.
				for i in range(n-1):
					if DotVecs(c.co-tp.co,tp.GetEdgeNormal())<0 or c.link==tp.link:
						ok=1
						break
					tp=tp.next
				if not ok:
					if DotVecs(c.co-tp.co,CrossVec2D(self.co-tp.co))<0 or c.link==tp.link:
						ok=1
				if not ok:
					#print "Cross Fail at ",c.co
					return self.next.nFill(n,tolerance,stop)
				prev=c
				c=c.next
			#print "no cross OK ------------------"

			#Test again that edge doesn't break out of loop
			if self!=lastpt.next:
				if (not Over180(-1.0*prev.GetVector(),self.GetVector(),prev.GetEdgeNormal(),self.GetEdgeNormal())) and DotVecs(lastpt.co-self.co,prev.GetEdgeNormal())<tolerance:
					#print "Breaks out of original poly (test 2)",lastpt.co,self.co,prev.next.co,prev.co
					return self.next.nFill(n,tolerance,stop)


			#make facet
			c=self
			facet=[]
			for i in range(n):
				#print "\t",c.co
				facet.append(c.link)
				prev=c
				c=c.next
			#do next
			self.next=prev

			res=self.Fill()
			if res:
				res.append(facet)
				return res
			#for fault tolerance, something is returned if we got this far
			return [facet]
				
		else:
			#Apparently there is nothing else in the loop, so we're done
			c=self
			facet=[]
			for i in range(n):
				facet.append(c.link)
				c=c.next
			return [facet]
			
	def Fill(self):
		#Maybe this would be better if it kept track of the loop size...
		#only or two vx? Then done.
		if self.next==self or self.next.next==self:
			#print self,"end"
			return []
		#only three faces?
		if self.next.next.next==self:
			#print self,"end3"
			return [[self.link,self.next.link,self.next.next.link]]
		#try to fill using Quads
		ret=self.nFill(4,0.01)
		if ret:
			#print self,"quads.."
			return ret
		#else try tris
		ret=self.nFill(3,0.01)
		if ret:
			#print self,"tris.."
			return ret 
		ret=self.nFill(4,0.0)
		if ret:
			#print self,"quads.. hightolerance"
			return ret
		#else try tris
		#print self,"high tolerance tris?"
		return self.nFill(3,0.0)
		
		
	def display(self,stop):
		print "\t",self.co,self.link.tco
		if self==stop:
			return
		self.next.display(stop)

class BSPElement2D:
	def __init__(self,pt,normal):
		self.pt=pt
		self.normal=normal
		self.inside=0
		self.outside=0
	def TestPoint(self,pt):
		return DotVecs(pt-self.pt,self.normal)
	def AddEdge(self,ptA,ptB,normal):
		testA=self.TestPoint(ptA)
		testB=self.TestPoint(ptB)
		if testA<0 or testB<0 or (testA==0 and testB==0):
			if self.outside:
				self.outside.AddEdge(ptA,ptB,normal)
			else:
				self.outside=BSPElement2D(ptA,normal)
		if testA>0 or testB>0  or (testA==0 and testB==0):
			if self.inside:
				self.inside.AddEdge(ptA,ptB,normal)
			else:
				self.inside=BSPElement2D(ptA,normal)
	def TestPointRecusive(self,pt):
		if self.TestPoint(pt)>0:
			if self.inside:
				return self.inside.TestPointRecusive(pt)
			else:
				return True
		else:
			if self.outside:
				return self.outside.TestPointRecusive(pt)
			else:
				return False

class VertexLoop:
	def __init__(self):
		self.loopBase=self.loopLast=0
		self.enclosed=[]
		self.BSPBase=0
		
	#Adding points to the loop
	def AddPoint(self,co,link):
		if self.loopLast:
			ne=LoopElement(co,link)
			self.loopLast.SetNext(ne)
			self.loopLast=ne
			self.bb=ExpandBounds2D(self.bb,co)
		else:
			self.loopBase=self.loopLast=LoopElement(co,link)
			self.bb=BuildBounds2D(co)
	def DonePoints(self):
		self.loopLast.SetNext(self.loopBase)
		self.loopBase=self.loopBase.next
		
	#Are the verts ordered in the right direction for an outer loop?
	def IsWrongDirection(self):
		
		#Find the leftmost point on the object
		leftmost=self.loopBase
		lp=leftmost.GetNextCo()
		cb=self.loopBase.GetNext()
		while cb!=self.loopBase:
			if cb.GetNextCo()[0]<lp[0]:
				leftmost=cb
				lp=cb.GetNextCo()
			cb=cb.GetNext()
		
		#Now compare it with a point even further left -
		#if the normals indicate that it is "inside", then the loop is backwards
		leftmostnext=leftmost.GetNext()
		refVect=Vector(self.bb[0][0]-1.0,self.bb[0][1])-leftmostnext.GetCo()
		if Over180((-1)*leftmost.GetVector(),leftmostnext.GetVector(),leftmost.GetEdgeNormal(),leftmostnext.GetEdgeNormal()):
			if DotVecs(refVect,leftmost.GetEdgeNormal())>0 or DotVecs(refVect,leftmostnext.GetEdgeNormal())>0:
				return True
			else:
				return False
		else:
			if DotVecs(refVect,leftmost.GetEdgeNormal())>0 and DotVecs(refVect,leftmostnext.GetEdgeNormal())>0:
				return True
			else:
				return False
			
	
	def Reverse(self):
		self.loopBase.Reverse(self.loopBase,0)

		
	#BSPs for determining loop heirarchy
	def BuildBSP(self):
		current=stop=self.loopBase
		self.BSPBase=BSPElement2D(current.GetCo(),current.GetEdgeNormal())
		while (current.GetNext()!=stop):
			current=current.GetNext()
			self.BSPBase.AddEdge(current.GetCo(),current.GetNextCo(),current.GetEdgeNormal())
	def TestLoopInLoop(self,other):
		if not TestBounds2D(self.bb,other.bb):
			return False
		if not self.BSPBase:
			self.BuildBSP()
		return self.BSPBase.TestPointRecusive(other.loopBase.GetCo())
	
	def TestEdgeIntersect(self,pA,pB):
		bb=BuildBounds2D(pA)
		bb=ExpandBounds2D(bb,pB)
		self.loopBase.GetNext().TestEdgeIntersect(pA,pB,bb,self.loopBase)
	
	def JoinLoop(self,other):
		#Reasonably big distance - bigger than anithing inside the shape
		dist=(self.bb[1][0]-self.bb[0][0])*(self.bb[1][1]-self.bb[0][1])
		pair=[]
		
		Ostop=Ocurrent=self.loopBase.GetNext()
		prevVect=self.loopBase.GetVector()*(-1)
		prevNorm=self.loopBase.GetEdgeNormal()
		Istop=other.loopBase
		while (1):
			
			Icurrent=Istop
			while (1):
				Icurrent=Icurrent.GetNext()
				
				vect=Icurrent.GetCo()-Ocurrent.GetCo()
				evect=Ocurrent.GetVector()
				enorm=Ocurrent.GetEdgeNormal()
				o180=Over180(evect,prevVect,enorm,prevNorm)

				#distance is less than current winner?
				d=DotVecs(vect,vect)
				#print "Distance: ",d
				if d<dist:
					#connecting edge would stay inside the outer loop?
					if (o180 and (DotVecs(vect,enorm)>0 or DotVecs(vect,prevNorm))) or ((not o180) and (DotVecs(vect,enorm)>0 and DotVecs(vect,prevNorm))):
						#And not intersect with anything?
						if not self.TestEdgeIntersect(Icurrent.GetCo(),Ocurrent.GetCo()):
							#And not with any other loops:
							ok=1
							for lp in self.enclosed:
								if lp.TestEdgeIntersect(Icurrent.GetCo(),Ocurrent.GetCo()):
									ok=0
									break
							if ok:
								#print "Winner: ",d
								dist=d
								pair=[Icurrent,Ocurrent]
							
				
				if Icurrent==Istop:
					break
			
			prevVect=evect*(-1)
			prevNorm=enorm
			Ocurrent=Ocurrent.GetNext()
			if Ocurrent==Ostop:
				break
	
		#success?
		if pair:
			#keep the edges going in the right direction...
			Icurrent=pair[0]
			Ocurrent=pair[1]
			#print "I/O",Icurrent.GetCo(),Ocurrent.GetCo()
			other.Reverse()
			#Connect the two loops
			Icopy=Icurrent.copy()
			Ocopy=Ocurrent.copy()
			Ocurrent.SetNext(Icopy)
			Icurrent.SetNext(Ocopy)
			return True
		print "ERROR: Could not connect loops."
		return False
	
	
	#Recursive filling!
	def Fill(self):
		facets=[]
		#print "----Fill Loops: ",self,facets
		#If there are loops encompased by this loop, they must be done first
		if (self.enclosed):
			#If there are any loops inside this loop, fill then first
			for eeloop in self.enclosed[0].enclosed:
				#print "Fill enc"
				facets.extend(eeloop.Fill())
					
			#Get on with the connecting...
			self.JoinLoop(self.enclosed[0])
			#print "Fill(): Connect:",self.loopBase.co,self.enclosed[0].loopBase.co
			self.enclosed.remove(self.enclosed[0])
			facets.extend(self.Fill())
			return facets
		#self.Display()
		res=self.loopBase.Fill()
		if res:
			facets.extend(res)
		return facets

#Figure out loop directions and such, building up loop objects
#Input: loops - an array of arrays of objects with a ".co" attribute
#       norm - Face normal
#Ouput: array of arrays or arrays [EdgeNormal,Coordinate,LinkToOriginal]
def BuildLoops(loops,norm):
	#Re-orient points so we can safely use just X&Y coordinates
	# (note: if Normal *dot* <0,0,1>  is >>0, this isn't really necissary)
	Basis=BuildBasis(norm,loops[0][1].tco-loops[0][0].tco)
	
	NewLoops=[]
	for loop in loops:
		
		#Build Loop Content
		NL=VertexLoop()
		for pt in loop:
			NL.AddPoint(Make2D(Basis,pt.tco),pt)
		NL.DonePoints()
		
		#NL.Reverse()
		
		#Fix direction if need be
		if NL.IsWrongDirection():
			print "loop at ",NL.loopBase.co," is backwards"
			NL.Reverse()

		NewLoops.append(NL)
		
	return NewLoops

#Make loops into a heirarchy of loops
def BuildHierarchy(loops):
	for loopA in loops:
		for loopB in loops:
			if loopA!=loopB:
				if loopA.TestLoopInLoop(loopB):
					loopA.enclosed.append(loopB)
					loopA.enclosed=BuildHierarchy(loopA.enclosed)
					loops.remove(loopB)
					return BuildHierarchy(loops)
	return loops
				

def BuildFacets(loops,normal):
	#Build up loops in the structure we need, and make sure
	#they are going the right direction
	newloops=BuildLoops(loops,normal)
	
	#if their are multiple loops, make a loop heirarchy
	if len(newloops)>1:
		newloops=BuildHierarchy(newloops)
	
	facets=[]
	for loop in newloops:
		r=loop.Fill()
		print r
		if r:
			facets.extend(r)
	return facets


def FillLoops(loops,oldFace,newMesh,flip,selected):
	if loops==[]:
		return
	fl=BuildFacets(loops,oldFace.normal*flip)
	for facet in fl:
		AddFace(oldFace,newMesh,facet,selected)

#Add a face to an nmesh
#Select verts as applicable
def AddFace(oldFace,newMesh,vertexList,selected):
	vxl=[]
	for v in vertexList:
		if v.newCounterpart==0:
			v.newCounterpart=NMesh.Vert(v.tco[0],v.tco[1],v.tco[2])
			newMesh.verts.append(v.newCounterpart)
			v.used=1
		if selected:
			v.newCounterpart.sel=1
		vxl.append(v.newCounterpart)
	newFace=NMesh.Face(vxl)
	if oldFace:
		if oldFace.counterPart:
			newFace.smooth=oldFace.counterPart.smooth
	newMesh.faces.append(newFace)

###################################################################
# Boolean Op objects                                              #
###################################################################

# 0 -> Unknown
# 1 -> Outside
# 2 -> Inside
# 3 -> Coplanar -> normals in same direction
# 4 -> Coplanar -> opposing normals
# 5 -> cut face
# 6 -> All done processing; ignore
		
class bVert:
	def __init__(self,coord,normal=0):
		global traversalCookie
		self.travCookie=traversalCookie
		coord.resize3D()
		self.tco=coord
		self.sharedvert=0
		if normal:
			self.normal=normal
		else:
			self.normal=Vector([0.0,0.0,0.0,1.0])
		self.e=[]
		self.f=[]
		self.newCounterpart=0
		self.redirect_to=0
		self.redirect_from=0

	def findSide(self,otherObject):
		print "v.fS(): \t\t\tlame vert side look up -",self,self.tco
		#Note: Maybe it would be faster to find the closeest face, and use its normal
		#      rather than doing an intersection test (_slightly_) less math
		ray=Vector([1,0,0])
		ifaces=0
		for f in otherObject.faces:
			if TestIntersectRayFace(self.tco,ray,f,1):
				ifaces+=1
		if (ifaces%2)==1:
			return 2
		else:
			return 1

	def GetRedirectTo(vA):
		if vA.redirect_to==0:
			return vA
		return vA.redirect_to.GetRedirectTo()
	def GetRedirectFrom(vA):
		if vA.redirect_from==0:
			return vA
		return vA.redirect_from.GetRedirectFrom()
	
	def real(self):
		return self
			
class bEdge:
	def __init__(self,bv0,bv1,infinite=0):
		global traversalCookie
		self.travCookie=traversalCookie
		i=0
		self.children=[]
		self.cutpoint=0
		self.f=[]
		self.v=[bv0,bv1]
		self.bb=BuildBounds(bv0.tco,bv1.tco)
		self.vect=bv1.tco-bv0.tco
		self.len=self.vect.length
		self.dot=DotVecs(self.vect,self.vect)
		self.infinite=0
		self.state=0
		self.sidescomputed=0
		self.broken=0
		self.linkto=0
		bv0.e.append(self)
		bv1.e.append(self)
		print "@@@@Adding Edge... ",self,bv0.tco,bv1.tco

	def setSides(self,state):
		#Maybe this should make sure it only operates on original edges
		if state==0:
			print "e.sS(): BUG: Tried to set side state to 0"
			return
		if self.sidescomputed:
			return
		print "e.sS()",state
		for f in self.f:
			f[0]=state
		self.sidescomputed=1

	def calcSides(self,force=0):
		#Maybe this should only be called onces per SuperMakeLoop?
		print "e.cS()",self
		if self.sidescomputed:
			return
		print "e.cS(): Compute Side: ",self
		if len(self.f)==2:
			if self.f[0][2]==self.f[1][2]:
				if self.f[0][0]:
					self.f[1][0]=self.f[0][0]
					self.sidescomputed=1
					return
				elif self.f[1][0]:
					self.f[0][0]=self.f[1][0]
					self.sidescomputed=1
					return
		for fA in self.f:
			if fA[0]==0 or fA[0]==5:
				closest=-1.1
				fv=fA[1]
				for fB in self.f:
					if fB[2].parent!=fA[2].parent:
						d=DotVecs(fA[1],fB[1])
						if d>closest:
							closest=d
							cf=fB
				if closest>-1.0:
					#print "e.cS(): ..."
					if closest>0.9999:
						print "e.cS(): COPLANAR FACES!!!!!"
						if DotVecs(fA[2].normal,cf[2].normal)>0:
							fA[0]=cf[0]=3
						else:
							fA[0]=cf[0]=4
					else:
						d=DotVecs(fA[1],cf[2].normal)
						#print "!!!!!!!!!!!this face:",fA,"other:",cf
						if (d>-0.0001 and d < 0.0001):
							print "e.cS(): Eeh??? BUG? Coplanar, but not?"
						elif d>0:
							fA[0]=1
						else:
							fA[0]=2
						#Can't judge other face's in/out automatically (probablems with "T" type setup, where one side has two coplanar faces)
						#d=DotVecs(cf[1],fA[2].normal)
						#if (d>-0.0001 and d < 0.0001):
						#	print "Eeh??? BUG? Coplanar, but not?"
						#elif d>0:
						#	cf[0]=1
						#else:
						#		cf[0]=2
				else:
					print "e.cS(): Need to determine side recursively..."
					if not force:
						return 0
		print "e.cS(): Computed Sides!",self.f		
		self.sidescomputed=1

				
	def testState(self,state,face):
		for f in self.f:
			#print "\tState: ",self,f[0],self.state,f[1],face
			if f[0]==state and f[1]==face:
				return 1
		return 0
		
			
class bFace:
	def __init__(self,nmFace,parent,virtual=0):
		global traversalCookie
		self.travCookie=traversalCookie
		
		self.parent=parent
		self.state=0
		self.e=[]
		self.v=[]
		self.extraVerts=[]
		self.allVerts=[]
		self.allEdges=[]
		self.center=Vector([0.0,0.0,0.0])
		self.altered=0
		self.built=0
		
		print "\tAdding Face..."
		self.counterPart=nmFace


		if virtual:
			self.normal=0
		else:
			if parent.matrix:
				self.normal=VecMultMat(Vector(list(nmFace.no)+[1]),parent.norm)
				self.normal.resize3D()
			else:
				self.normal=Vector(nmFace.no)
			self.normal.normalize()
			
			#First the verts
			for v in nmFace.v:
				bv=parent.verts[v.index]
				bv.f.append(self) #maybe not needed...
				self.v.append(bv)
				self.allVerts.append(bv)
				self.center+=bv.tco
			self.center/=len(self.v)
			
			self.bb=BuildBounds(self.v[0].tco,self.v[1].tco)
			self.bb=ExpandBounds(self.bb,self.v[2].tco)
			
			self.e.append(MakeEdge(self.v[0],self.v[1],self))
			self.e.append(MakeEdge(self.v[1],self.v[2],self))
			
			if len(self.v)==3:
				self.e.append(MakeEdge(self.v[2],self.v[0],self))
			else:
				self.bb=ExpandBounds(self.bb,self.v[3].tco)
				self.e.append(MakeEdge(self.v[2],self.v[3],self))
				self.e.append(MakeEdge(self.v[3],self.v[0],self))
				

	def buildEdges(self):
		#print "f.bE():\tBuilding new edges for ",self,self.altered
		if self.altered:
			c=len(self.allVerts)
			print "f.bE():\tBuilding new edges for ",self," Verts:",c
			for A in range(c):
				for B in range(A+1,c):
					vA=self.allVerts[A]
					vB=self.allVerts[B]
					go=0 #Should an edge be made?
					mo=1 #Should this data be merged with an existing edge
					for fA in vA.f:
						if fA.parent!=self.parent:
							for fB in vB.f:
								if fB.parent!=self.parent:
									if fA==fB:
										d=DotVecs(fA.normal,self.normal)
										if (d>-0.9999 and d<0.9999):
											mo=0
										go=1
					if go:
						print "f.bE():\t\tNew Edge: ",A,vA.tco,B,vB.tco
						ne=MakeEdge(vA,vB,self,1,mo)
						#if ne:
						#	ne.calcSides()
	
						
	
	def constructFace(self,otherObject,destMesh,inside,outside,cosame,coopposing,flip=1):
		#Should detect unaltered/simple quad/tri faces faster
		global enorm
		print "f.cf():\tConstructing face ",self
		if self.built:
			return
		#for v in self.allVerts:
		#	print "f.cf():\t\t",v.tco
		loops=[[],[],[],[],[]]
		for e in self.allEdges:
			e.calcSides()
		for e in self.allEdges:
			sides=[]
			for f in e.f:
				if f[2]==self:
					sides.append([f[0],f[1]])
			print "f.cf():\t\tE: ",e,e.v[0].tco,e.v[1].tco,sides
		for e in self.allEdges:
			print "f.cf():\t\tEdge: ",e,e.v[0].tco,e.v[1].tco
			if (abs(e.v[0].tco[0]-0.75)<0.1 and abs(e.v[1].tco[0]-1.0)<0.1 and abs(e.v[0].tco[1]+1.00)<0.1):
				print "#############################################"
				print e.f
			if (abs(e.v[1].tco[0]-0.75)<0.1 and abs(e.v[0].tco[0]-1.0)<0.1 and abs(e.v[0].tco[1]+1.00)<0.1):
				print "#############################################"
				print e.f
			#e.calcSides()
			if enorm==1:
				AddFace(self,destMesh,e.v,0)
					
			for f in e.f:
				#if 0:
				#draw edge normal
				pA=MidPoint(e.v[0].tco,e.v[1].tco)
				pB=pA+f[1]*0.05
				if enorm==1:
					AddFace(self,destMesh,[NewVx(pA),NewVx(pB)],0)
					if len(e.children)>0 or e.linkto or e	.broken:
						AddFace(self,destMesh,[NewVx(e.v[0].tco+f[1]*0.02),NewVx(e.v[1].tco+f[1]*0.02)],0)

				if f[2]==self and f[3]==0 and (f[0]>0 or f[0]<5):
					#NOTE: Maybe by using the edge normal to determine to use v[0] or v[1]
					#could be used to automatic normals
					l = SuperMakeLoop(e,e.v[0],e.f.index(f),e)
					if len(l):
						if l[0]>4:
							l[0]=4
						print "f.cf():loop of type ",l[0]
						if l[0]==0:
							#Manually figure out side...
							l[0]=e.v[0].findSide(otherObject)
							print "f.cf(): Lame Lookup side: ",l[0]
						#There should problably be a better way of negotiating the side while traversing the edges...
						for le in l[1]:
							if le[0].broken==0:
								le[0].setSides(l[0])
						loops[l[0]].append(l[1])
						for q in l[1]:
							print "f.cf(): \t",q[1].tco
							if enorm==2:
								AddFace(self,destMesh,q[0].v,1)
								for ff in q[0].f:
									if ff[2]==self:
										pA=MidPoint(q[0].v[0].tco,q[0].v[1].tco)
										pB=pA+ff[1]*0.05
										AddFace(self,destMesh,[NewVx(pA),NewVx(pB)],0)
					else:
						print "f.cf(): LOOP FAILED"
		if len(loops)==0:
			print "f.cf(): COULD NO FILL FACE",self
			return
		#NOTE: This can probably be done a little more efficiently now
		#that edge normals are precomputed and stored
		vxloops=[[],[],[],[],[],[]]
		for l in range (1,5):
			print "f.cf(): Loops of type ",l
			for lp in loops[l]:
				print "f.cf(): \tNL"
				nloop=[]
				for li in lp:
					nloop.append(li[1])
					print "f.cf(): \t\t",li[1].tco
				vxloops[l].append(nloop)
		if outside:
			#print "Fill outside"
			FillLoops(vxloops[1],self,destMesh,flip,0)
		if inside:
			#print "Fill inside"
			FillLoops(vxloops[2],self,destMesh,flip,0)
		if cosame:
			#print "Fill cosame"
			FillLoops(vxloops[3],self,destMesh,flip,0)
		if coopposing:
			#print "Fill coop"
			FillLoops(vxloops[4],self,destMesh,flip,0)
		self.built=1
		#Go do other faces next
		for e in self.e:
			if e.broken==0:
				for f in e.f:
					if f[2].parent==self.parent and f[2].built==0:
						print "f.cf(): continuing to neighbor ",f
						f[2].constructFace(otherObject,destMesh,inside,outside,cosame,coopposing,flip)
						

class bMesh:
	def __init__(self,blenderObject,matrix):
		print "Adding mesh from: %s"%blenderObject.getName()
		
		#setup variables
		self.verts=[]
		self.faces=[]
		self.edges=[]
		self.infinite=0
		
		m=NMesh.GetRawFromObject(blenderObject.getName())
		self.nmesh=m
		
		#process verts
		self.matrix=matrix
		if (matrix):
			self.rot=rotmatrix=matrix.rotationPart()
			self.norm=CopyMat(matrix)
			self.norm.invert()
			self.norm.transpose()
			for v in m.verts:
				self.verts.append(bVert(VecMultMat(Vector(list(v.co)+[1]),matrix),VecMultMat(Vector(list(v.no)+[1]),self.norm)))
		else:
			for v in m.verts:
				self.verts.append(bVert(Vector(list(v.co)),Vector(list(v.no))))
					
		#process faces
		for f in m.faces:
			if len(f.v)>2:
				self.faces.append(bFace(f,self))
			else:
				#stuff related to infinite faces
				1
		#if there are just edges & not faces, make this face object infinite

	def constructFaces(self,otherObject,destMesh,inside,outside,cosame,coopposing,flip=1):
		print "m.cF(): Constructing Faces for ",self
		t=time.time()
		for f in self.faces:
			if f.altered:
				f.constructFace(destMesh,otherObject,inside,outside,cosame,coopposing,flip)
		#return
		#take care of any loose ends in a 2nd pass
		for f in self.faces:
			if not f.altered:
				f.constructFace(destMesh,otherObject,inside,outside,cosame,coopposing,flip)
		print "Face fill time: ",time.time()-t

def MakeEdge(v0,v1,face,inout=0,mergeonly=0):
	
	if v0==v1:
		print "MakeEdge(): BUG! Attempt to make edge pointing to self."
		return []
	cve=v1.tco-v0.tco
	if cve.length<0.0001:
		print "MakeEdge(): BUG! Overlapping verticies",v0,v1,v0.tco,v1.tco
		return []
	
	#Check if it already exists...
	ENorm=CrossVecs(face.normal,v1.tco-v0.tco)
	#ENorm=CrossVecs(v1.tco-v0.tco,face.normal)
	ENorm.normalize()
	
	for e in v0.e:
		if e.v[0]==v1 or e.v[1]==v1:
			for f in e.f:
				if f[2]==face:
					print "MakeEdge(): BUG??? Tried to add edge to face that already had it..."
					return e
			e.f.append([0,ENorm,face,0])
			if inout:
				e.f.append([0,ENorm*(-1),face,0])
			face.allEdges.append(e)
			return e
	#if not, make a new edge
	if not mergeonly:
		e=bEdge(v0,v1)
		e.f.append([0,ENorm,face,0])
		if inout:
			e.f.append([0,ENorm*(-1),face,0])
		face.parent.edges.append(e)
		face.allEdges.append(e)
		return e
	return []
	

###################################################################
# Merging Tools to keep things clean                              #
###################################################################

def MergeVert(vA,vB):
	if vA==vB:
		return vA
	v=vA.tco-vB.tco
	#print "###################################################################"
	#print "MERGEVERTICIES: ",v.lenth,vA.tco,vB.tco
	for f in vB.f:
		if not vA in f.allVerts:
			f.allVerts.remove(vB)
			f.allVerts.append(vA)
			vA.f.append(f)
		if vB in f.v:
			f.v[f.v.index(vB)]=vA
		if vB in f.parent.verts:
			f.parent.verts[f.parent.verts.index(vB)]=vA
	for e in vB.e:
		if vB in e.v:
			e.v[e.v.index(vB)]=vA
			vA.e.append(e)
	for f in vA.f:
		f.altered=1
	#test for merging edges
	for eA in vA.e:
		for eB in vA.e:
			if eA!=eB:
				if eA.v[0]==vA:
					rA=eA.v[1]
				else:
					rA=eA.v[0]
				if eB.v[0]==vA:
					rB=eB.v[1]
				else:
					rB=eB.v[0]
				if rA==rB:
					MergeEdge(eA,eB)
				dv=rA.tco-rB.tco
				if dv.length<0.0001:
					MergeVert(rA,rB)
				#if (eA.v[0]==eB.v[0] and eA.v[1]==eB.v[1]) or (eA.v[0]==eB.v[1] and eA.v[1]==eB.v[0]):
				#	#print "EDGE OVERLAP",eA,eB
	return vA


def MergeEdge(eA,eB):
	if eA==eB:
		return eA
	for v in eB.v:
		v.e.remove(eB)
	for f in eB.f:
		while eB in f[2].allEdges:
			f[2].allEdges.remove(eB)
		if not eA in f[2].allEdges:
			f[2].allEdges.append(eA)
		if eB in f[2].e:
			f[2].e[f[2].e.index(eB)]=eA
		while eB in f[2].parent.edges:
			f[2].parent.edges[f[2].parent.edges.index(eB)]=eA
		eA.f.append([0,f[1],f[2],0])
	#eA.broken=1 #mark edge as altered
	eB.broken=1	
	eB.linkto=eA
	eB.f=[] # this shouldn't be necissary, but is for one of the test cases
	        # it is probably indicative of a bug somewhere else
	return eA
	

###################################################################
# Functions for Applying Intersection                             #
###################################################################

def xptdraw(nm):
	global xpt,drawxpt
	if drawxpt:
		for pt in xpt:
			if pt[0]==0:
				p=pt[1]
				AddFace(0,nm,[NewVx(p+Vector(0.0,0.01,0.0)),NewVx(p-Vector(0.0,0.01,0.0))],0)
				AddFace(0,nm,[NewVx(p+Vector(0.01,0.0,0.0)),NewVx(p-Vector(0.01,0.0,0.0))],0)
				AddFace(0,nm,[NewVx(p+Vector(0.0,0.0,0.01)),NewVx(p-Vector(0.0,0.0,0.01))],0)
			if pt[0]==1:
				p=pt[1]
				AddFace(0,nm,[NewVx(p+Vector(0.007,0.007,0.0)),NewVx(p-Vector(0.007,0.007,0.0))],0)
				AddFace(0,nm,[NewVx(p+Vector(0.007,0.0,0.007)),NewVx(p-Vector(0.007,0.0,0.007))],0)
				AddFace(0,nm,[NewVx(p+Vector(0.0,0.007,0.007)),NewVx(p-Vector(0.0,0.007,0.007))],0)

#new vertex for a given coordinate
def NewVx(co):
	return bVert(co)

def LinkVx(vA,vB):
	MergeVert(vA,vB)
	return 1

def CutEdge(e,co):
	#note: If an edge is cut, we can probably immediately mark inside vs outside
	#      Problems only arise when edges overlap, which happens with vx in edge / vx in vx situations
	if e.linkto:
		print "CutEdge(): BUGBUGBUG!!!! Tried to cut inactive edge"
		return CutEdge(e.linkto,co)
	cv=co-e.v[0].tco
	#note: must scale by length so fuzzy matching is kept sane
	#note: maybe scaled point could be passed along, so there is no need
	#      to do the dot product over & over again
	cutat=cv.length
	if (cutat<-0.0001 or cutat>e.len+0.0001):
		print "CutEdge(): BUG! Tried to cut edge out of bounds"
	
	if e.children:
		cd=cutat-e.cutPoint
		if cd<0.0001 and cd>-0.0001:
			#same point
			return e.children[0].v[1]
		elif cd>0:
			return CutEdge(e.children[1],co)
		else:
			return CutEdge(e.children[0],co)
	
	#test for cut point being at one of the two endpoints	
	if cutat<0.0001:
		return e.v[0]
	elif cutat>(e.len-0.0001):
		return e.v[1]
	
	e.broken=1
	#Edge actually needs to be partitioned
	print "CutEdge(): Partitioning Edge",e,"at",co,"(",cutat,"/",e.len,")"
	nv=NewVx(co)
	e.v[0].e.remove(e)
	e.v[1].e.remove(e)
	e.cutPoint=cutat
	A=bEdge(e.v[0],nv)
	B=bEdge(nv,e.v[1])
	e.children.append(A)
	e.children.append(B)
	for f in e.f:
		nv.f.append(f[2])
		f[2].allEdges.remove(e)
		f[2].allEdges.append(A)
		f[2].allEdges.append(B)
		B.f.append([0,f[1],f[2],0])
		A.f.append([0,f[1],f[2],0])
		f[2].allVerts.append(nv)
		f[2].altered=1

	return nv

def VxInFace(f,v):
	if not f in v.f:
		v.f.append(f)
	if not v in f.allVerts:
		f.allVerts.append(v)
	for vf in v.f:
		vf.altered=1
	return 1

def VxInEdge(e,v):
	vp=CutEdge(e,v.tco)
	LinkVx(vp,v)
	return 1


###################################################################
# The Magic - the Intersection detector & processor               #
###################################################################

fxf=0
exf=0
vxf=0
exe=0
vxe=0
vxv=0

#sfxf=0
sexf=0
svxf=0
sexe=0
svxe=0
svxv=0

fexf=0
ffxf=0
fexe=0

def IntersectVertexVertex(vA,vB):
	global vxv,svxv
	vxv+=1
	#d0=vA.tco[0]-vB.tco[0]
	#d1=vA.tco[1]-vB.tco[1]
	#d2=vA.tco[2]-vB.tco[2]
	#d=d0*d0+d1*d1+d2*d2
	d=vA.tco-vB.tco
	if d.length<0.0001:
		#Verticies are effectively the same
		print "-VxV() Vertex overlap at ",vA.tco,vA,vB
		svxv+=1
		LinkVx(vA,vB)
		for f in vA.f:
			f.altered=1
		return 1
	return 0

def IntersectVertexEdge(v,e):
	global vxe,svxe
	vxe+=1
	i=0
	#print "VxE",v.tco,e.v[0].tco,e.v[1].tco
	#Actually, maybe this loop should be _after_ TestVertInEdge3D
	for v2 in e.v:
		i+=IntersectVertexVertex(v,v2)
	if (i==0) and (TestVertInEdge3D(v,e)):
		print "-VxE() Vertex in edge at ",v.tco
		svxe+=1
		nv=CutEdge(e,v.tco)
		LinkVx(nv,v)
		i=1
	return i

def IntersectEdgeEdge(eA,eB):
	global exe,fexe,sexe,xpt
	exe+=1
	if TestBounds(eA.bb,eB.bb):
		fexe+=1
		i=0
		A=eA.v[1].tco-eA.v[0].tco
		B=eB.v[1].tco-eB.v[0].tco
		#used to be ,0,0 at end - edit if intersecting edges don't work
		#Note: Tolerance should be smaller
		r=FindIntersectLineLine3D(eA.v[0].tco,eA.v[1].tco,eB.v[0].tco,eB.v[1].tco,1,1,1,0.000001)
		if len(r)>2:
			p=eA.v[0].tco+r[0]*A
			#test for vertex in edge
			for v in eA.v:
				i+=IntersectVertexEdge(v,eB)
			for v in eB.v:
				i+=IntersectVertexEdge(v,eA)
			if (i==0):
				print "-ExE(): Edge intersects edge at ",p
				print "\tQualification: ",r
				print "\nEdges: ",eA.v[0].tco,eA.v[1].tco,eB.v[0].tco,eB.v[1].tco
				xpt.append([1,p])
				sexe+=1
				nvA=CutEdge(eA,p)
				nvB=CutEdge(eB,p)
				LinkVx(nvA,nvB)
				i=1
		else:
			A.normalize()
			B.normalize()
			dv=DotVecs(A,B)
			if (dv>0.999 or dv <-0.999): #parallel?
				print "ExE(): Edges parallel",A,B
				for v in eA.v:
					i+=IntersectVertexEdge(v,eB)
				for v in eB.v:
					i+=IntersectVertexEdge(v,eA)
		return i
	return 0

def IntersectVertexFace(v,f):
	#Note vertex may have already been added to
	#the face due to an intersection with an edge
	global vxf,svxf,xpt
	vxf+=1
	#print "VxF",v.tco,TestPointInFace(v.tco,f)
	if TestPointInFace(v.tco,f) and TestPointInPlane(v.tco,f.normal,f.v[0].tco):
		i=0
		for v2 in f.allVerts:
			i+=IntersectVertexVertex(v,v2)
		if (i==0):
			#Add point to face
			print "-VxF(): Point in face at ",v.tco,f,
			svxf+=1
			VxInFace(f,v)
			xpt.append([1,v.tco])
		return 1
	return 0

def IntersectEdgeFace(e,f):
	global exf,fexf,sexf,xpt
	exf+=1
	if TestBounds(e.bb,f.bb):
		fexf+=1
		#Intersect?
		p=FindIntersectEdgeFace(e,f)
		if p:
			#print "FxE",p
			i=0
			for eB in f.e:
				i+=IntersectEdgeEdge(e,eB)
			for v in e.v:
				i+=IntersectVertexFace(v,f)
			#iff vert doesn't lie in the object's face, we need to cut the edge
			if (i==0):
				#cut edge
				#note: at this point we know inside vs outside iff there isn's something screwy
				#      further down the edge
				print "-ExF(): Face cuts Edge at ",p
				sexf+=1
				nv = CutEdge(e,p)
				VxInFace(f,nv)
				xpt.append([0,nv.tco])
				i=1
			return i
				
		else:
			#Parallel?
			i=0
			dv=DotVecs(e.v[0].tco-e.v[1].tco,f.normal)
			if (dv<0.00001 and dv>-0.00001):
				#print "Edge Parllel to face"
				for eB in f.e:
					i+=IntersectEdgeEdge(e,eB)
				for v in e.v:
					i+=IntersectVertexFace(v,f)
			return i
	return 0

def FullIntersect(meshA,meshB):
	#is this FacexFace method faster, or all edges vs all edges in two passes?
	global traversalCookie,fxf,ffxf
	traversalCookie+=1
	travBase = traversalCookie

	print meshA,meshB
	
	#Iterate Face vs Face
	for fA in meshA.faces:
		traversalCookie+=1
		fATrav = traversalCookie
		for fB in meshB.faces:
			traversalCookie+=1
			fBTrav = traversalCookie
			fxf+=1
			if TestBounds(fA.bb,fB.bb):
				ffxf+=1
				#print "Possibly intersecting faces: ",fA,fB
				#fA.e vs fB -> can eliminate wth travcookie, test ExF, E.vxF, ExE, ExV, VxV
				print "eA vs fB"
				for eA in fA.e:
					if 1 or eA.travCookie<fATrav: # travcookie works for exf, but not vxf and such
						IntersectEdgeFace(eA,fB)
						eA.travCookie=traversalCookie
				#fB.e vs fA -> must test all, only ExF?
				print "eB vs fA"
				for eB in fB.e:
					if 1 or eB.travCookie<fATrav:
						IntersectEdgeFace(eB,fA)
						eB.travCookie=traversalCookie
				#Note: after FxF, maybe it would make sense to build new edges immediately to avoid looking through long lists for matches...
	for fA in meshA.faces:
		fA.buildEdges()
	for fB in meshB.faces:
		fB.buildEdges()

###################################################################
# Base functions direction the operation                          #
###################################################################


def CookieCutter(cookie,cutter):
	x=Draw.PupMenu("Faces to keep%t|Keep all|Inside only|Outside Only")
	if x<1:
		return

	MatA=Matrix(*map(list,cookie.getMatrix()))
	MatB=Matrix(*map(list,cutter.getMatrix()))
	#MatB.invert()
	#MatA=MatA*MatB
	#MatB.identity()

	bCookie=bMesh(cookie,MatA)
	bCutter=bMesh(cutter,MatB)

	#IntersectMeshes(bCutter,bCookie,1)
	#print "###############################################################"
	#IntersectMeshes(bCookie,bCutter,0)
	#IntersectMeshes(bCookie,bCutter,1)
	#print "###############################################################"
	#IntersectMeshes(bCutter,bCookie,0)
	FullIntersect(bCookie,bCutter)
	
	nm=NMesh.New()
	
	if x==1:
		bCookie.constructFaces(nm,bCutter,1,1,1,1)
	elif x==2:
		bCookie.constructFaces(nm,bCutter,1,0,1,0)
	elif x==3:
		bCookie.constructFaces(nm,bCutter,0,1,0,1)

	nm.mode=bCookie.nmesh.mode
	NMesh.PutRaw(nm,"foo"	)
	return

def Difference(cookie,cutter):
	MatA=Matrix(*map(list,cookie.getMatrix()))
	MatB=Matrix(*map(list,cutter.getMatrix()))
	#MatB.invert()
	#MatA=MatA*MatB
	#MatB.identity()

	bCookie=bMesh(cookie,MatA)
	bCutter=bMesh(cutter,MatB)

	#IntersectMeshes(bCookie,bCutter,1)
	#IntersectMeshes(bCutter,bCookie,0)
	FullIntersect(bCookie,bCutter)
	
	nm=NMesh.New()
	bCutter.constructFaces(nm,bCookie,1,0,0,1,-1)
	bCookie.constructFaces(nm,bCutter,0,1,0,0,1)

	nm.mode=bCookie.nmesh.mode
	NMesh.PutRaw(nm,"foo")
	return

def Intersection(cookie,cutter):
	MatA=Matrix(*map(list,cookie.getMatrix()))
	MatB=Matrix(*map(list,cutter.getMatrix()))
	#MatB.invert()
	#MatA=MatA*MatB
	#MatB.identity()

	bCookie=bMesh(cookie,MatA)
	bCutter=bMesh(cutter,MatB)

	#IntersectMeshes(bCookie,bCutter,1)
	#IntersectMeshes(bCutter,bCookie,0)
	FullIntersect(bCookie,bCutter)
	
	nm=NMesh.New()
	bCutter.constructFaces(nm,bCookie,1,0,1,0)
	bCookie.constructFaces(nm,bCutter,1,0,0,0)
	#bCookie.constructFaces(nm,bCutter,0,1,0,0)
	xptdraw(nm)
	
	nm.mode=bCookie.nmesh.mode
	NMesh.PutRaw(nm,"foo")
	
	return

def Union(cookie,cutter):
	MatA=Matrix(*map(list,cookie.getMatrix()))
	MatB=Matrix(*map(list,cutter.getMatrix()))
	#MatB.invert()
	#MatA=MatA*MatB
	#MatB.identity()

	bCookie=bMesh(cookie,MatA)
	bCutter=bMesh(cutter,MatB)

	#IntersectMeshes(bCookie,bCutter,1)
	#IntersectMeshes(bCutter,bCookie,0)
	FullIntersect(bCookie,bCutter)
	
	nm=NMesh.New()
	bCutter.constructFaces(nm,bCookie,0,1,1,0)
	bCookie.constructFaces(nm,bCutter,0,1,0,0)

	nm.mode=bCookie.nmesh.mode
	NMesh.PutRaw(nm,"foo")
	
	return

def DoBool():
	#Check the sanity of the selection
	objsel=Object.GetSelected()
	meshes=0
	for ob in objsel:
		if ob.getType() == 'Mesh':
			meshes += 1
	if meshes!=2 :
		Draw.PupMenu("ERROR: Exactly two mesh objects must be selected")
		return
	x=Draw.PupMenu("MegaBool%t|1: Intersect|2: Union|3: Difference|4: Cookie Cutter")
	#print x
	if x==1:
		Intersection(objsel[1],objsel[0])
	elif x==2:
		Union(objsel[1],objsel[0])
	elif x==3:
		Difference(objsel[1],objsel[0])
	elif x==4:
		CookieCutter(objsel[1],objsel[0])

DoBool()
#oA=Object.Get("Cookie")
#oB=Object.Get("Cutter")
#oB=Object.GetSelected()[0]
#print oA
#print oB
#Intersection(Object.Get("Cookie"),Object.Get("Cutter"))
#Difference(oA,oB)
#Union(oA,oB)
#CookieCutter(Object.Get("Cookie"),Object.Get("Cutter"))
#CookieCutter(oB,oA)

print "Done. Tests:\n\tFxF",fxf,"\t",ffxf,"\n\tExF",exf,"\t",fexf,"\t",sexf
print "\tVxF",vxf,"\t--\t",svxf,"\n\tExE",exe,"\t",fexe,"\t",sexe,"\n\tVxE",vxe,"\t--\t",svxe,"\n\tVxV",vxv,"\t--\t",svxv
			
if len(messages)>0 and 0:
	txt=Blender.Text.New("Boolean Messages")
	txt.clear()
	txt.write(messages)


#Note to self:
# Re-orient the messhes when doing a 2D cut so
#two-axis bounding box test can be done
#difference op should maybe just flip the normals of the B object and do an intersection
#
#For f-gons, don't cut with internal edges
#when making face do whole f-gon at once (so outer edge & inner loops)