#!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-2004  Rewrite of alot - completely changed how inside   <br>
                vs outside determination works                   <br>
                                                                 
"""

__author__ = "Theodore K. Schundler"
__version__ = "r10 8-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 *

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 ...

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)
						#note: the nsl>0.0001 thing is for dealing with edges where ang is 180 degrees (-1)
						if (ang>0 or nsl>0.0001) and (abs(dQ-dW)<0.001) and ang>bestAngle:
							#print "ML():\t",nsum,dA,dB
							#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                                              #
###################################################################


#This may be a little faster if conversted to a 2D function using the desired
#normal as the Z axis (doing a matrix transformation)

#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)

def JoinConcentric(vertList,enclosed,CCWlist,allList,normal,f,nm):
	#Preffer parallel edges (smaller dot product absolute value)
	#something is wrong here... should be both close _and_ parallel
	#not just close...
	bestOuter=0
	bestInner=0
	parallelness=10000000
	Ocount=len(vertList)
	for On in range(Ocount):
		I0=vertList[(On+Ocount-1)%Ocount].tco-vertList[On].tco;
		I1=vertList[(On+1)%Ocount].tco-vertList[On].tco;
		for vl in enclosed:
			Icount=len(vl)
			for In in range(Icount):
				vect=vertList[On].tco-vl[In].tco;
				d=vect.length;
				par=parallelness
				O0=vl[(In+Icount-1)%Icount].tco-vl[In].tco;
				O1=vl[(In+1)%Icount].tco-vl[In].tco;
				pt=(abs(DotVecs(I0,O0))+0.01)*d
				if (pt<par):
					par=pt
				pt=(abs(DotVecs(I0,O1))+0.01)*d
				if (pt<par):
					par=pt
				pt=(abs(DotVecs(I1,O0))+0.01)*d
				if (pt<par):
					par=pt
				pt=(abs(DotVecs(I1,O1))+0.01)*d
				if (pt<par):
					par=pt
				#is this worth considering?
				vc=vertList[On].tco-vl[In].tco
				par=par*DotVecs(vc,vc)
				if (par<parallelness):
					#does this not intersect w/ any existing edges?
					if TestEdgeIntersectLoop3D(vertList[On],vl[In],vertList)==0:
						ok=1
						for vl2 in CCWlist:
							if TestEdgeIntersectLoop3D(vertList[On],vl[In],vl2):
								#print On," to ",On," failed. ",vl,vl2
								ok=0
								break
						if ok:
							bestOuter=On
							bestInner=(vl,In)
							parallelness=par
	if parallelness<10000000:
		#if we have a winner, add its points to the real loop
		#this means one edge is in there twice
		#(as it should be since it is used by two faces)
		newVL=[]
		#print "Connect ",bestOuter," to ",bestInner[1]
		for n in range(len(vertList)):
			v=vertList[n]
			if n==bestOuter:
				newVL.append(v)
				print v.tco,vertList.index(v)
				vl=bestInner[0]
				for i in range(bestInner[1],len(vl)):
					newVL.append(vl[i])
					#print vl[i].tco,":",i
				for i in range(bestInner[1]+1):
					newVL.append(vl[i])
					#print vl[i].tco,":",i
				newVL.append(v)
				print "Join Concentric():",v.tco,vertList.index(v)
			else:
				print "JoinConcentric():",v.tco,vertList.index(v)
				newVL.append(v)
		#remove the added points from the lists of loops
		#that might be inside this loop
		CCWlist.remove(bestInner[0])
		if FillProcessedLoop(newVL,CCWlist,allList,normal,f,nm):
			return 1
		else:
			print "JoinConcentric(): attempt to connect loops failed... (uh-oh... bug?)"
	else:
		print "JoinConcentric(): couldn't find path to connect concentric loops"
	return 0

#Note: for testing concentric, maybe just find closest (don't need to find real length, just len^2 via dot product) edge and use its normal?

#fill in terms of n-sided convex polygons
def nFillLoop(sides,vertList,CCWlist,normal,oldFace,newMesh,selected,tolerance):
		#try all ways to fill this loop using points on the edges
		#there are len(vertList) possible faces we need to test
		count=len(vertList);
		for i in range(count):
			convex=1
			#test convexity
			loop=[]
			loopReal=[]
			for j in range(sides):
				d=DotVecs(CrossVecs(vertList[(i+(j+1)%sides)%count].tco-vertList[(i+j)%count].tco,vertList[(i+(j+sides-1)%sides)%count].tco-vertList[(i+j)%count].tco),normal)
				#print d,0.000001
				if d<tolerance:
					convex=0
					break
				loop.append((i+j)%count)
				loopReal.append(vertList[(i+j)%count])
			#if its convex, then we can investigate further
			#print "CX:",convex,loop,vertList
			if convex:
				#test no verts in face
				ok=1
				for p in vertList:
					nocheck=0
					for lp in loopReal: #don't test the verts that make the new face
						if p==lp:
							nocheck=1

					if nocheck==0 and TestPointInConvexLoop(p.tco,loopReal,normal):
						ok=0
						break
				#if the new edge doesn't cut into existing edges
				#(determined with the TestPointInConvexLoop above)
				if ok:
					ok=1
					#test verts in inner loop
					enc=[]
					for vl in CCWlist:
						for lp in vl:
							if TestPointInConvexLoop(lp.tco,loopReal,normal):
								try:
									enc.index(vl)
								except:
									enc.append(vl)
								ok=0
					if ok:
						newVL=[]
						n=loop[sides-1]
						while n!=loop[1]:
							newVL.append(vertList[n])
							n+=1
							n%=count
						if FillProcessedLoop(newVL,CCWlist,normal,oldFace,newMesh,selected):
							#print "Face! ",[vertList[i],vertList[(i+1)%count],vertList[(i+2)%count],vertList[(i+3)%count]]
							AddFace(oldFace,newMesh,loopReal,selected)
							return 1
						else:
							1
							#print "failed..."
							#for vx in vertList:
							#	print vx.tco
							#print "--------------------------"
					else:
						#There is apparently a heirarchy of loops here
						#Connect this loop to a revelvant inner loop
						return JoinConcentric(vertList,enc,CCWlist,normal,oldFace,newMesh,selected)
			#else:
			#	print "Not Convex..."
		#print "Fail fill size",sides
		return 0


def FillProcessedLoop(vertList,CCWlist,normal,oldFace,newMesh,selected):
	#Four or more points, then we try to fill with quads
	if len(vertList)>3:
		if nFillLoop(4,vertList,CCWlist,normal,oldFace,newMesh,selected,0.00001):
			return 1
	#Three points or quad fill failed, we try to fill with tris
	if len(vertList)>2:
		if nFillLoop(3,vertList,CCWlist,normal,oldFace,newMesh,selected,0.00001):
			return 1
		else:
			#if it failed, don't be as picky about what is convex.
			if len(vertList)>3:
				if nFillLoop(4,vertList,CCWlist,normal,oldFace,newMesh,selected,0):
					return 1
				#Three points or quad fill failed, we try to fill with tris
			if nFillLoop(3,vertList,CCWlist,normal,oldFace,newMesh,selected,0):
					return 1
			else:
					return 0
	#if there are less than three points we're done
	else:
		return 1


def FillPreprocess(VertexLists,normal):
	RightWay=[]
	WrongWay=[]

	for ThisList in VertexLists:
		lastTest=ThisList[len(ThisList)-1]
		FlipList=0
		Direction=0
		#this runs the test multiple times. Really just two times where you get the same result in
		#a row should be sufficient, or once in "fast mode"
		for testpt in ThisList:
			#Have the ray start midway on the segment because it may get confused on the edges
			midpt=MidPoint(lastTest.tco,testpt.tco)
			#The ray is perpendicular to the edge and the normal
			ray=CrossVecs(testpt.tco-lastTest.tco,normal)
			#print "B:",ThisList[0].tco,ThisList[1].tco,"->",midpt,ray

			isect_counter_all=0
			isect_counter_others=0

			for OtherList in VertexLists:
				prevVert=OtherList[len(OtherList)-1]
				for Vert in OtherList:
					#print Vert.tco
					if not ((prevVert==lastTest)): #Don't even try to test for intersection with self
						#print "TX: ",midpt,midpt+ray,prevVert.tco,Vert.tco,testpt.tco,lastTest.tco
						isect=FindIntersectLineLine3D(midpt,midpt+ray,prevVert.tco,Vert.tco,1,0,0)
						#if len(isect):
						#	print "\tX:",isect[0],isect[1],midpt+ray*isect[0],isect[2]
						if (len(isect)>0 and isect[0]>0):
							isect_counter_all+=1
							if (not OtherList==ThisList):
								isect_counter_others+=1
					prevVert=Vert

			#if the ray intersects an odd number of edges, it is pointing in the wrong direction
			#if it is pointing in the wrong direction, the verticies aren't ordered right
			if (isect_counter_all&1==1):
				#print "Flipped List"
				FlipList+=1
				#ThisList.reverse()
			else:
				FlipList-=1

			#if it intersects an even number of edges of objects other than it self
			#then it is an outer loop in a set of concentric loops (i.e. the outer edge of a doughnut)
			if (isect_counter_others&1==0):
				Direction+=1
				#RightWay.append(ThisList); 
			#otherwise it is an inner loop (the edges that make up the hole of a doughnut)
			#and the verts are ordered in the other direction
			else:
				Direction-=1
				#WrongWay.append(ThisList);
				
			lastTest=testpt
			
		#apply direction...
		#print "Flipped: ",FlipList,"Direction:",Direction 
		if FlipList>0:
			ThisList.reverse()
		if Direction>0:
			RightWay.append(ThisList)
		else:
			WrongWay.append(ThisList)
		lastTest=testpt
		
	#print "Right Way:",len(RightWay),"WrongWay:",len(WrongWay)
	return [RightWay,WrongWay]


def FillLoops(loops,oldFace,newMesh,flip,selected):

	if loops==[]:
		return
		
	#print "SRART FILLL------"		
		
	#if the face needs to be flipped (i.e. when doing a difference operation)
	normal=oldFace.normal*flip

	#First establish loop directions
	FixedLists=FillPreprocess(loops,normal)	

	#Now make two lists - clockwise and counterclockwise
	CWlist=FixedLists[0]
	CCWlist=FixedLists[1]

	#Time to fill!
	for l in CWlist:
		i=FillProcessedLoop(l,CCWlist,normal,oldFace,newMesh,selected)
		if i==0:
			print "FillLoops():\t\t\tFill Failed..."
			for v in l:
				print "FillLoops():\t\t\t\t",v.tco,v


###################################################################
# 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(nmFace.no),parent.rot)
			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,4):
			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:
			FillLoops(vxloops[1],self,destMesh,flip,0)
		if inside:
			FillLoops(vxloops[2],self,destMesh,flip,0)
		if cosame:
			FillLoops(vxloops[3],self,destMesh,flip,0)
		if coopposing:
			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()
			for v in m.verts:
				self.verts.append(bVert(VecMultMat(Vector(list(v.co)+[1]),matrix),VecMultMat(Vector(list(v.no)),rotmatrix)))
		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
		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)

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.vect
		B=eB.vect
		#used to be ,0,0 at end - edit if intersecting edges don't work
		r=FindIntersectLineLine3D(eA.v[0].tco,eA.v[1].tco,eB.v[0].tco,eB.v[1].tco,1,1)
		if len(r)==3:
			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
				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"
				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(Object.Get("Cookie"),Object.Get("Cutter"))
#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)