#!BPY
# -*- coding: latin-1 -*-
"""
Name: 'Isobaren'
Blender: 248
Group: 'Mesh'
Tooltip: 'subdivide parallel to axis'
"""

__author__ = ["Mathias Weitz"]
__url__ = ("blender")
__version__ = "0.0.2 17/3/08"

__bpydoc__ = """\
The script adds lines in a mesh that are parallel to the x, y, z-levels are

Es gibt zwei Parameter für den Script. There are two parameters for the script.

1.) Die Ebene, an denen die Isobaren ausgerichtet werden sollen 

1) The level at which the isobars should be targeted

2.) Abstand der Isobaren 

2) spacing of the isobars

Die Isobaren können nicht entlang der Achsen verschoben werden, die Ausrichtung ist immer am Ursprung. 

The isobars can not be moved along the axis, the orientation is always at the origin. 

Man kann ja sein Mesh verschieben und auch drehen, daß es passt - das ist einfacher und übersichtlicher als weitere Parameter. I

t may well be moving mesh and rotate it fits - it is simpler and clearer than other parameters.

...
"""
# -------------------------------------------------------------------------- 
# Isobaren by Mathias Weitz
# -------------------------------------------------------------------------- 
# ***** BEGIN GPL LICENSE BLOCK ***** 
# 
# This program is free software; you can redistribute it and/or 
# modify it under the terms of the GNU General Public License 
# as published by the Free Software Foundation; either version 2 
# of the License, or (at your option) any later version. 
# 
# This program is distributed in the hope that it will be useful, 
# but WITHOUT ANY WARRANTY; without even the implied warranty of 
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
# GNU General Public License for more details. 
# 
# You should have received a copy of the GNU General Public License 
# along with this program; if not, write to the Free Software Foundation, 
# Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA. 
# 
# ***** END GPL LICENCE BLOCK ***** 
# -------------------------------------------------------------------------- 

import Blender
from Blender import *
from Blender import NMesh
from Blender.BGL import *
from Blender.Draw import *

import math
from math import *

# Events
EVENT_NOEVENT	= 1
EVENT_START	  = 2
EVENT_EXIT	   = 3
EVENT_STARTTEST  = 4
EVENT_DIST	   = 5
EVENT_X		  = 6
EVENT_Y		  = 7
EVENT_Z		  = 8
EVENT_MODE	   = 9
EVENT_MARKFACES  = 10
EVENT_REDUCEFACES  = 11

# Inits
V_DIST = 1.0
V_X = 1
V_Y = 1
V_Z = 1
V_MODE = 1

######################################################
# GUI drawing
######################################################
def draw():
	global EVENT_NOEVENT,EVENT_START,EVENT_EXIT,EVENT_MARKFACES,EVENT_REDUCEFACES 
	global EVENT_DIST, V_DIST, B_DIST
	global EVENT_X, V_X, B_X
	global EVENT_Y, V_Y, B_Y
	global EVENT_Z, V_Z, B_Z
	global EVENT_MODE, V_MODE, B_MODE

	#global EVENT_STARTTEST
	
	
	########## Titles
	glClear(GL_COLOR_BUFFER_BIT)
	glColor3f(0,0,0)
	textPy = 270
	glRasterPos2d(10,textPy)
	Text ("Isobaren-Divide v0.0.2")
	textPy -= 15
	glRasterPos2d(15,textPy)
	Text ("- Choose several Faces and press 'start isobaren'")
	textPy -= 15
	glRasterPos2d(15,textPy)
	Text ("- Reducing Faces only fill Faces, if there is")
	textPy -= 15
	glRasterPos2d(35,textPy)
	Text ("just one outline with four or less Verts")
	textPy -= 15
	glRasterPos2d(35,textPy)
	Text ("otherwise just the outline ist created")
	textPy -= 15
	#
	glColor3f(0.1, 0.5, 0.1)
	glRecti(5, 40, 240, 95)
	
	glColor3f(0.1, 0.1, 0.5)
	glRecti(5, 100, 240, 195)

	B_X = Toggle("X", EVENT_X, 10, 170, 75, 18, V_X)
	B_Y = Toggle("Y", EVENT_Y, 85, 170, 75, 18, V_Y)
	B_Z = Toggle("Z", EVENT_Z, 160, 170, 75, 18, V_Z)

	B_DIST = Slider("Edge: ",EVENT_DIST , 10, 150, 225, 18, V_DIST, 0.05, 5.0, 1)

	B_MODE = Menu ("Resolve torsion %t|Torsion Quads to Triangles %x1|ZYX %x2|YXZ %x3",EVENT_MODE, 10, 130, 225, 18 ,V_MODE)		 
	Button("Start Isobaren",EVENT_START , 10, 110, 225, 18)

	Button("Start select same faces", EVENT_MARKFACES , 10, 70, 225, 18)
	Button("Start reduce faces", EVENT_REDUCEFACES , 10, 50, 225, 18)

	Button("Exit",EVENT_EXIT , 10, 10, 225, 18)

	# Button("StartTest",EVENT_STARTTEST , 10, 30, 120, 18)
		

def event(evt, val):
	if (evt == QKEY and not val): 
		Exit()

def bevent(evt):
	global EVENT_NOEVENT,EVENT_START,EVENT_EXIT,EVENT_MARKFACES,EVENT_REDUCEFACES 
	#global EVENT_STARTTEST
	global EVENT_DIST, V_DIST, B_DIST
	global EVENT_X, V_X, B_X
	global EVENT_Y, V_Y, B_Y
	global EVENT_Z, V_Z, B_Z
	global EVENT_MODE, V_MODE, B_MODE
	
	if evt == EVENT_EXIT:
		Exit()	  
	elif evt == EVENT_X:
		V_X = B_X.val
	elif evt == EVENT_Y:
		V_Y = B_Y.val
	elif evt == EVENT_Z:
		V_Z = B_Z.val
	elif evt == EVENT_DIST:
		V_DIST = B_DIST.val
	elif evt == EVENT_MODE:
		V_MODE = B_MODE.val
	elif evt == EVENT_REDUCEFACES:
		reduceFaces()
	elif evt == EVENT_MARKFACES:
		markFaces()
		Blender.Redraw()
	elif evt == EVENT_START:
		startDraw(V_DIST,[V_X,V_Y,V_Z],V_MODE)
		Blender.Redraw()
	
#################################
# Die Routinen fuer Points
#################################

def projPointA(p1,p2,p3,pt,dir,mode):
	"berechnet die senkrechten bzw durch dir(ection) bestimmte Projektion vom KoordinatenArray p \
	 auf der Ebene, die durch p2-p1 und p3-p1 aufgespannt wird \
	 Rueckgabe ist ein Array[xyz] der Koordinaten auf der Ebene \
	 und True, falls der Punkt in der Ebene liegt. \
	 mode = 1 wird auf die Kontrolle bei dreieckigen Flaechen verzichtet \
	 damit bekommt man viereckige Flaechen besser im Griff, aber es gibt auch mehr Fehler"   
	erg = [pt[0], pt[1], pt[2], True,0]
	abw = 0.001
	d1x = p2[0] - p1[0]
	d1y = p2[1] - p1[1]
	d1z = p2[2] - p1[2]
	d2x = p3[0] - p1[0]
	d2y = p3[1] - p1[1]
	d2z = p3[2] - p1[2]
	ddx = pt[0] - p1[0]
	ddy = pt[1] - p1[1]
	ddz = pt[2] - p1[2]
	if dir[0] == 0 and dir[1] == 0 and dir[2] == 0:
		ec = crossProd([d1x,d1y,d1z], [d2x,d2y,d2z])
	else:
		ec = dir
	e =  norm(ec)
	m = matrix3x3([d1x,d1y,d1z], [d2x,d2y,d2z],e)
	m.t()   
	detm = m.det()
	if 0.000001 < abs(detm):
		so =m.solve([ddx, ddy, ddz])
		b = True
		c1x = p2[0]*so[0] + p1[0]*(1-so[0])
		c1y = p2[1]*so[0] + p1[1]*(1-so[0])
		c1z = p2[2]*so[0] + p1[2]*(1-so[0])
		c2x = p3[0]*so[1] + p1[0]*(1-so[1])
		c2y = p3[1]*so[1] + p1[1]*(1-so[1])
		c2z = p3[2]*so[1] + p1[2]*(1-so[1])
		# Test, ob Punkt in der Ebene liegt
		if so[0] < -abw or (1+abw) < so[0] or so[1] < -abw or (1+abw) < so[1]:
			b = False
		if mode == 0 and (1+abw) < abs(so[0]) + abs(so[1]):
			b = False 
		erg = [c1x+c2x-p1[0], c1y+c2y-p1[1], c1z+c2z-p1[2], b, so[2]]
	
	return erg  


def projPoint(p1,p2,p3,pt,dir,mode):
	return projPointA(pointToArray(p1), pointToArray(p2), pointToArray(p3), \
		pointToArray(pt), dir, mode)
	 

def projPointEA(faceP, p, dir=[0,0,0]):
	"berechnet die senkrechten Projektion vom KoordinatenArray p \
	 auf der Ebene. Bei einer Ebene mit 4 Punkten wird \
	 der Punkt und der Abstand gemittelt."
	erg = [0,0,0,False,0]
	c = len(faceP.v)
	#print c,
	if (3 < c):
		# bei mehr als 3 Ecken wird das Ergebnis gemittelt   
		sumx = 0;
		sumy = 0;
		sumz = 0;
		sumdir = 0
		hit = False
		rc = 0
		for i in range(c):
			im = (i-1)%c
			ip = (i+1)%c
			pjt = projPointA( \
				pointToArray(faceP.v[i]), \
				pointToArray(faceP.v[im]), \
				pointToArray(faceP.v[ip]), \
				p, dir ,0)
			if pjt[3]:
				rc += 1
				sumx += pjt[0]
				sumy += pjt[1]
				sumz += pjt[2]
				sumdir += pjt[4]
			hit |= pjt[3]
			
		if 0 < rc:
			erg = [sumx / rc, sumy / rc, sumz / rc, hit, sumdir / rc]
	else:
		pjt = projPointA( \
			pointToArray(faceP.v[0]), \
			pointToArray(faceP.v[1]), \
			pointToArray(faceP.v[2]), \
			p, dir ,0)
		erg = [pjt[0], pjt[1], pjt[2], pjt[3], pjt[4]]
	return erg 

def projPointE(faceP,p, dir=[0,0,0]):
	return projPointEA(faceP, pointToArray(p), dir)

def edgeDistA (a1s,a1e,a2s,a2e,d = [0,0,0]):
	"wie edgeDistP, nur es werden 3-Array statt Verts uebergeben"
	erg = [0,0,0,False]
	# e ist das normierte Kreuzprodukt
	
	d1 = vecSub(a1e, a1s)
	d1neg = vecSub(a1s, a1e)
	d2 = vecSub(a2e, a2s)
	d2neg = vecSub(a2s, a2e)
	dd = vecSub(a1s, a2s)   
	ddneg = vecSub(a2s, a1s)
	de = d
	if d == [0,0,0]:
		de = crossProd(d1,d2)
	e = norm(de)
	m = matrix3x3(d1neg,d2,e)
	if abs(m.det()) == 0:
		if d == [0,0,0]:
			# beide Geraden laufen parallel
			# es gibt eventuell viele Loesungen
			# TODO
			pass
		else:
			# es Geraden laufen aus Sicht von d parallel
			pass
	else:
		m.t()
		so = m.solve(dd)
		erg[0] = so[0]
		erg[1] = so[1]
		erg[2] = so[2]
		if 0 <= so[0] and so[0] <= 1 and 0 <= so[1] and so[1] <= 1:
			erg[3] = True
	return erg  

def edgeDistP (p1s,p1e,p2s,p2e,d = [0,0,0]):
	"Kuerzeste Distanz zweier Geraden \
	  das Ergebniss ist ein Array mit \
		0. m1 = Verhaeltniss auf der ersten Geraden \
		1. m2 = Verhaeltniss auf der zweiten Geraden \
		2. dist = Distanz der beiden Geraden \
		3. b = True, wenn die kuerzeste Distanz innerhalb der Kanten liegt"
	return edgeDistA (pointToArray(p1s),pointToArray(p1e),pointToArray(p2s),pointToArray(p2e),d)

def edgeDistQ (p1s,p1e,p2s,p2e):
	"Suche naechsten Punkt zur Geraden \
	 p1e-p1s und p2e-p2s \
	 Ergebniss ist ein Array mit\
	   1. der Punkt auf der Geraden p2e-p2s \
	   2. True, falls sich der Punkt innerhalb von p1e-p1s befindet"
	erg = [0,False]
	z =  edgeDistP (p1s,p1e,p2s,p2e)
	mo = z[0]
	mb = False
	if 0 <= mo and mo <= 1:
		mb = True
	m = z[1]
	if m < 0:
		m = 0
	if 1 < m:
		m = 1
	erg[0] = propEdge(p2s, p2e, m)
	erg[1] = mb
	return erg 

def findClosestFace(me, co, dir):
	"Sucht das Face, das in der Richtung dir dem \
	Koordinaten in co am naechsten ist \
	Rueckgabe [face, coProj, dist]"
	T_Dir = 0
	coProj = [0,0,0]
	minAbw = 0.00001
	bestFace = 0
	besth = 10000   
	for ff in me.faces:
		pjt1 = projPointEA(ff,co, dir)
		if pjt1[3]:											 
			d = dist([co[0] - pjt1[0], co[1] - pjt1[1], co[2] - pjt1[2]])
			if minAbw < d and d < besth and (pjt1[4] <= 0 or (pjt1[4] != 0 and T_Dir == 1)):
				besth = d
				bestFace = ff
				coProj[0] = pjt1[0]
				coProj[1] = pjt1[1]
				coProj[2] = pjt1[2]
	if bestFace != 0:
		pass
		#pc1 = NMesh.Vert(coProj[0], coProj[1], coProj[2])
		#pc2 = NMesh.Vert(co[0], co[1], co[2])
		#me.verts.append(pc1)
		#me.verts.append(pc2)
		#e = me.addEdge(pc1,pc2)
		#po[i][1] = bestFace
		#po[i][2] = pc1
	else:
		print "found no bestFace: ", besth
	return [bestFace, coProj, besth]

def findAllEdges(me, v1, v2, dir):
	"Gibt eine Liste aller Edges zurueck \
	 die Listenelemente enthaelt folgende Infos \
	 0. ProjEdge \
	 1. Koordinaten auf der Ausgangsgeraden \
	 2. Distanz"
	crossEdges = []
	mminAbw = 0.0001
	for ee in me.edges:
		vv1 = pointToArray(ee.v1)
		vv2 = pointToArray(ee.v2)
		if not arrayEqual(v1, vv1) and not arrayEqual(v1, vv2) and not arrayEqual(v2, vv1) and not arrayEqual(v2, vv2): 
			er = edgeDistA(v1, v2, vv1, vv2,dir)
			if er[3] and 0 < er[2]:
				# Test, ob Kante verdeckt wird durch ein Face
				bestPoint1 = ratioPoint (v1, v2, er[0])
				erf = findClosestFace(me, bestPoint1, dir)
				if er[2] < erf[2] + mminAbw:
					crossEdges.append(ee)
	
	erg = [[0,v1,0],[0,v2,0]]
	b = True
	sec = 0
	while b:
		b = False
		newerg = []
		newerg.append(erg[0])
		for i in range(len(erg)-1):
			best = 100000
			targetEdge = 0
			bestRatio1 = 0
			bestPoint1 = [0,0,0]
			bestPoint2 = [0,0,0]
			for ee in crossEdges:
				vv1 = pointToArray(ee.v1)
				vv2 = pointToArray(ee.v2)
				er = edgeDistA(erg[i][1], erg[i+1][1], vv1, vv2,dir)
				if er[3]:	   
					if er[2] < best and mminAbw < er[0] and er[0] < 1 - mminAbw:
						best = er[2]
						targetEdge = ee
						bestRatio1 = er[0]
						bestPoint1 = ratioPoint (erg[i][1],erg[i+1][1],er[0])   
						bestPoint2 = ratioPoint (vv1,vv2,er[1])
		
			if best < 100000:   
				#print sec, bestRatio1, best	
				#pc1 = NMesh.Vert(bestPoint1[0], bestPoint1[1], bestPoint1[2])
				#pc2 = NMesh.Vert(bestPoint2[0], bestPoint2[1], bestPoint2[2])
				#me.verts.append(pc1)
				#me.verts.append(pc2)
				#e = me.addEdge(pc1,pc2)
				newerg.append([targetEdge, bestPoint1, best, bestPoint2])
				sec += 1
				if sec < 3500:
					b = True
			newerg.append(erg[i+1])
		erg = newerg
	return erg
		
				
def propEdge(p1, p2, m):
	"Gibt die Verhaeltnisskoordinaten \
	 zwischen den beiden Punkten zurueck"
	return [p2.co.x*m + p1.co.x*(1-m), p2.co.y*m + p1.co.y*(1-m), p2.co.z*m + p1.co.z*(1-m)]

def findFaces(me, ar):
	"die Faces, das die angegebenen Punkte enthaelt, \
	rueckgabe ist ein Array mit den faces \
	die Komplexitaet ist |f|*|ar|*|aveEdges| also etwa Anzahl der Faces * 16"
	global T_Ign
	erg = []
	for f in me.faces:
		if T_Ign != 1 or f.hide != 1:
			b = True
			for j in ar:
				found = False
				for k in f.v:
					if j == k:
						found = True
				b &= found
			if b:
				erg.append(f)
	return erg

def findFaces2(me, ar):
	"wie oben, etwas einfacher und ohne T_Ign"
	erg = []
	for f in me.faces:
		b = True
		for j in ar:
			found = False
			for k in f.v:
				if j == k:
					found = True
			b &= found
		if b:
			erg.append(f)
	return erg

def getConnected(me, v):
	"gibt eine Liste der verbundenen Punkte zurueck"
	erg = []
	for ee in me.edges:
		if ee.v1 == v:
			erg.append(ee.v2)
		if ee.v2 == v:
			erg.append(ee.v1)
	return erg

def getOpenConnected(me, v):
	"gibt eine Liste der verbundenen Punkte mit offener Kante zurueck"
	erg = []
	for ee in me.edges:
		b = False
		if ee.v1 == v:
			p = ee.v2
			b = True
		if ee.v2 == v:
			p = ee.v1
			b = True
		if b:
			e = findFaces2(me,[v,p])
			if len(e) < 2:
				erg.append(p)
	return erg

########################
# Such-Pfade
########################

class pfadElem:
	"Element"
	def __init__(self,vert,countc=0):
		self.v = vert
		self.count = countc
		self.pred = []
		# interne Verwaltung
		self.hasChanged = True  

	def append(self,node):
		"fuegt den pfadElem in die Nachbarliste \
		wird dadurch der Pfad verkurzt, dann setze count und gib True zurueck"
		erg = False
		b = True
		for vv in self.pred:
			if vv == node:
				b = False
			if node.count < self.count - 1:
				# der Knotenzaehler wurde veraendert
				self.count = node.count + 1
				erg = True
				self.hasChanged = True
		if b:
			self.pred.append(node)
		return erg
		
	def equalvert(self,vert):
		return self.v==vert
		
	def getRawPath(self):
		"erstellt vom Knoten zum Suchursprung einen Array mit den Verts"
		pass

	def show(self):
		print self.v.index , ":", self.count
	
def findOpenPath (me,v1,v2, wantedLength = 0):
	"sucht den kuerzesten Pfad von v1 nach v2 nur ueber offene Kanten"
	m = []
	erg = []
	m.append(pfadElem(v1))
	b = True
	# fuer Test, um sicher zu gehen, dass die Routine abbricht
	r = 50000
	while b and 0 < r:
		r = r -1 
		#print "************ next **********"
		#print "m = " , len(m)
		b = False
		#for v in m:
		#   v.show()
		#print "************ einzeltest **********" , r
		for v in m:
			#print dir(v)
			#print v.__dict__
			if v.hasChanged:
				v.hasChanged = False
				#v.show()
				w = getOpenConnected(me, v.v)
				#print "nachbarn", len(w)
				for vv in w:
					# Test ob Nachbarpunkt in Pfadmenge
					bb = True
					for a in m:
						#print "Test punkt bereits enthalten:" , vv.index , a.v.index,
						if a.equalvert(vv):
							bb = False
							#print "bereits enthalten",
							# die Counter updaten
							if v.v != v2:
								# ueber den Zielpunkt soll keine Suche hinausgehen
								b |= a.append(v)
							if a.v != v2:
								b |= v.append(a)
						#print
					if bb and 0 < r and v.v != v2:
						# ueber den Zielpunkt v2 soll keine Suche hinausgehen
						r = r - 1
						newa = pfadElem(vv,100000)
						b |= newa.append(v)
						v.append(newa)
						m.append(newa)
	
	# Ergebniss erstellen
	# start
	vgo = v2
	vg = 0
	for v in m:
		if v.equalvert(vgo):
			vg = v
			erg.append(vg.v)
	# wenn eine bestimmte Laenge gesucht
	if 0 < wantedLength:
		next = 0
		if vg != 0:
			for vn in vg.pred:
				#print "cmp" , vn.count, wantedLength
				if vn.count == wantedLength-1:
					next = vn
			if next != 0:
				vg = next
				erg.append(vg.v)
		
	while vg != 0:
		next = 0
		for vn in vg.pred:
			if vn.count < vg.count:
				next = vn
		vg = 0
		if next != 0:
			vg = next
			erg.append(vg.v)
			#vg.show()
		
	#print "Testlaenge", r
	#print 
	#print erg
	return erg  

def notIn(set1, set2):
	"Ergebniss ist ein Set von Elementen, die in Set1 aber nicht in Set2 sind"
	erg = []
	for v in set1:
		b = True
		for w in set2:
			if v == w:
				b = False
		if b:
			erg.append(v)
	return erg  

def crossProdP(p, p1, p2):
	"berechnet das Kreuzprodukt zwischen p1-p und p2-p"
	v1 = vecSub(pointToArray(p),pointToArray(p1))
	v2 = vecSub(pointToArray(p),pointToArray(p2))
	return crossProd(v1,v2)

def cosP(p, p1, p2):
	"berechnet den Cosinus p1-p und p2-p"
	v1 = vecSub(pointToArray(p),pointToArray(p1))
	v2 = vecSub(pointToArray(p),pointToArray(p2))
	return vecCos(v1,v2)

def queueIndexDist(i1,i2, plen):
	"die Distanz zweier Indezes in einer zyklischen Liste"
	min = abs(i1-i2)
	if i1 + plen - i2 < min:
		min = i1 + plen - i2
	if i2 + plen - i1 < min:
		min = i2 + plen - i1
	return min

def minA (li):
	erg = 0
	if 0 < len(li):
		erg = li[0]
	for w in li:
		if w < erg:
			erg = w
	return erg

def maxA (li):
	erg = 0
	if 0 < len(li):
		erg = li[0]
	for w in li:
		if erg < w:
			erg = w
	return erg

def arrayEqual(a1,a2):
	erg = False
	if len(a1) == len(a2):
		erg = True
		for i in range(len(a1)):
			if 0.0000001 < abs(a1[i] - a2[i]):
				erg = False
	return erg

def ratioPoint(a1,a2,r):
	return [a1[0]*(1-r)+a2[0]*r, a1[1]*(1-r)+a2[1]*r, a1[2]*(1-r)+a2[2]*r]


######################################
# mathematische Routinen
######################################

class matrix3x3:
	"Wird zum losen linearer Systeme der Grobe 3 benutzt. \
	Parameter sind 3 Array mit je 3 Elementen, \
	die wichtigste Routine ist solve(Point)"
	def __init__(self,p1,p2,p3):
		self.p1 = p1
		self.p2 = p2
		self.p3 = p3

	def det(self):
		"Determinante nach der sarrusschen Regel"
		return \
			self.p1[0] * (self.p2[1] * self.p3[2] - self.p2[2] * self.p3[1]) \
			- self.p2[0] * (self.p1[1] * self.p3[2] - self.p3[1] * self.p1[2]) \
			+ self.p3[0] * (self.p1[1] * self.p2[2] - self.p2[1] * self.p1[2])

	def t(self):
		h = self.p1[1]
		self.p1[1] = self.p2[0]
		self.p2[0] = h
		h = self.p1[2]
		self.p1[2] = self.p3[0]
		self.p3[0] = h
		h = self.p2[2]
		self.p2[2] = self.p3[1]
		self.p3[1] = h
		return self

	def write(self):
		print "[(" , self.p1[0] , ","  , self.p1[1] , "," , self.p1[2] , \
			") (" , self.p2[0] , ","  , self.p2[1] , "," , self.p2[2] , \
			") (" , self.p3[0] , ","  , self.p3[1] , "," , self.p3[2] ,  ")]"
		return self

	def solve(self,e):
		self.t()
		m = self.det()
		mx = matrix3x3(e, self.p2, self.p3)
		my = matrix3x3(self.p1, e, self.p3)
		mz = matrix3x3(self.p1, self.p2, e)
		erg = [0, 0, 0]
		if m != 0:
			erg = [mx.det()/m, my.det()/m, mz.det()/m]
		return erg

def printArray(v):
	print "[",
	for i in range(len(v)):
		print "%5.3f" % v[i],
	print "]",

def pointToArray(p1):
	return [p1.co.x, p1.co.y, p1.co.z]

def pointPrint(p1):
	return "[" + "%5.2f"%p1.co.x + "," + "%5.2f"%p1.co.y + "," + "%5.2f"%p1.co.z + "]"

def vecAdd(v1,v2):
	"subtraktion"
	return [v2[0]+v1[0], v2[1]+v1[1], v2[2]+v1[2]]

def vecSub(v1,v2):
	"subtraktion"
	return [v2[0]-v1[0], v2[1]-v1[1], v2[2]-v1[2]]
	
def scalarMul(v,a):
	"multiplikation mit einem Skalar"
	return [a*v[0], a*v[1], a*v[2]]

def propEdgeV(v1, v2, m):
	"Gibt die Verhaeltnisskoordinaten \
	 zwischen den beiden Punkten zurueck"
	return [v2[0]*m + v1[0]*(1-m), v2[1]*m + v1[1]*(1-m), v2[2]*m + v1[2]*(1-m)]

def dist(v1):
	"quadratische Norm"
	return math.sqrt(v1[0]*v1[0]+v1[1]*v1[1]+v1[2]*v1[2])

def norm(v1):
	"normalisierter Norm"
	d = dist(v1)
	if d == 0: d = 1
	return [v1[0]/d,v1[1]/d,v1[2]/d]

def scalarProd(v1,v2):
	"zwei 3-Arrays, es wird das Scalarprodukt berechnet"
	return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]

def crossProd(v1,v2):
	"zwei 3-Arrays, es wird das Kreuzprodukt berechnet"
	return [v1[1]*v2[2]-v1[2]*v2[1],v1[2]*v2[0]-v1[0]*v2[2],v1[0]*v2[1]-v1[1]*v2[0]]

def vecCos(v1,v2):
	"berechnet den Cosinus des Winkels zwischen zwei Vektoren"
	d1 = dist(v1)
	d2 = dist(v2)
	if d1 == 0:
		d1 = 1.0
	if d2 == 0:
		d2 = 1.0
	return scalarProd(v1,v2) / (d1 * d2)

def cosBetweenPoints(p0,p1,p2):
	a0 = pointToArray(p0)
	a1 = pointToArray(p1)
	a2 = pointToArray(p2)
	return vecCos(vecSub(a0,a1),vecSub(a0,a2)) 

def sinBetweenPoints(p0,p1,p2):
	a0 = pointToArray(p0)
	a1 = pointToArray(p1)
	a2 = pointToArray(p2)
	return crossProd(vecSub(a0,a1),vecSub(a0,a2)) 

def projPointV(k1,k2,pj):
	"Projektion von pj auf die durch k1,k2 aufgespannte Ebene \
	 Ergebniss ist der Punkt, die Verhaeltnisse, die Hoehe, \
	 und ein Flag, falls pj auf die Ebene projeziert wird"   
	erg = True
	
	ec = crossProd(k1, k2)
	e =  norm(ec)
	m = matrix3x3(k1, k2, e)
	m.t()
	so =m.solve(pj)
	
	cx = k1[0]*(1-so[0]) + k2[0]*(1-so[1])
	cy = k1[1]*(1-so[0]) + k2[1]*(1-so[1])
	cz = k1[2]*(1-so[0]) + k2[2]*(1-so[1])
	rp = [cx, cy, cz]
	# Test, ob Punkt in der Ebene liegt
	if so[0] < 0 or so[1] < 0:
		erg = False
		
	return [rp,so[0],so[1],so[2],erg] 


def edgeProj(startPunkt, m, endPunkt):
	"Gibt den Projektionspunkt von m auf die Gerade |startPunkt, endPunkt|"
	p = vecSub(startPunkt, endPunkt)
	v = vecSub(startPunkt, m)
	d = dist(p)
	if d == 0:
		d = 1.0
	return vecAdd(startPunkt, scalarMul(norm(p), scalarProd(p,v) / d))

def line(me, v1, v2):
	p1 = NMesh.Vert(v1[0], v1[1], v1[2])
	p2 = NMesh.Vert(v2[0], v2[1], v2[2])
			
	me.verts.append(p1)
	me.verts.append(p2)
	me.addEdge(p1, p2)
		
def lineDir(me, v1, d1):
	p1 = NMesh.Vert(v1[0], v1[1], v1[2])
	v2 = vecAdd(v1, d1)
	p2 = NMesh.Vert(v2[0], v2[1], v2[2])
			
	me.verts.append(p1)
	me.verts.append(p2)
	me.addEdge(p1, p2)
	
def polyCenter(poly):
	m = [0,0,0]
	for p in poly:
		pa = pointToArray(p)
		m = vecAdd(m, pa)   
		
	m = scalarMul(m,1.0 / len(poly))
	return m

def drawPoly(me, poly, mode = 0):
	m = polyCenter(poly)
	mm = newVert(me, m[0], m[1], m[2])
	
	if mode == 0:
		for p in poly:
			me.addEdge(mm, p)
	else:
		for pi in range(len(poly)):
			p1 = poly[pi]
			p2 = poly[(pi + 1) % len(poly)]
			pmm = scalarMul(vecAdd(p1,p2), 0.5)
			p = newVert(me, pmm[0], pmm[1], pmm[2])
			me.addEdge(mm, p)


def distAA(p1,p2):
	"Distanz zwischen zwei Punkten"
	return dist(vecSub(p1,p2))

def distPP(p1,p2):
	"Distanz zwischen zwei Punkten"
	return distAA(pointToArray(p1),pointToArray(p2))

##########################
# Punkte und Linien
##########################
	
def distPoint2Edge(p, p1, p2):
	return distAA2Edge(pointToArray(p), pointToArray(p1), pointToArray(p2))
	
def distAA2Edge(p, p1, p2):
	# ermittelt die minimale Distanz des Punktes p
	# zu der Linie p1, p2
	# Ergebniss
	# 0. Lotdistanz (skalar)
	# 1. Distanz zur Linie (skalar), falls das Lot nicht auf der Linie liegt wird die 
	#	Distanz zu einem Punkt gemessen.
	# 2. Lot liegt auf der Linie (boolean)
	# 3. Lotpunkt
	# 4. = 1 Lotpunkt liegt auf p1, = 2 Lotpunkt liegt auf p2, = 0 sonst
	
	v1 = vecSub(p1,p)
	v2 = vecSub(p2,p)
	d = vecSub(p1,p2)
	
	r1 = dist(v1)
	r2 = dist(v2)
	dd = dist(d)
	# ldist ist die Lotdistanz
	ldist = 0
	# adist ist die absolute Distanz zur Linie
	adist = 0
	# lot auf Linie
	blot = True
	mm = p
	atEnd = 0
	if 0 < dd:
		# m ist das Verhaeltniss entlang der markierten Geraden
		m = scalarProd(v1,d) / (dd*dd)
		if abs(m) < 0.00001:
			atEnd = 1 
		elif abs(m-1) < 0.00001:
			atEnd = 2 
		# mm ist der Lotpunkt
		mm = vecAdd(p1 , scalarMul(d, m))
		# mmb ist die Lotdistanz 
		ldist = dist(vecSub(p,mm))
		adist = ldist
		if (0 > m):
			# ausserhalb der Linie, naeher an p1
			adist = dist(v1)
			blot = False
		if (1 < m):
			# ausserhalb der Linie, naeher an p2
			adist = dist(v2)
			blot = False
		
		#print "Punkt ",
		#printArray(p)
		#print
		#print "Linie ",
		#printArray(p1)
		#printArray(p2)
		#print
		#print "Verhaeltniss auf der Linie ", m
		#print "Lotpunkt: ",
		#printArray(mm)
		#print
		#print "Lotdistanz: ", ldist
		#print "Absolute Distanz: ", adist
		#print
		
	return [ldist,adist,blot,mm, atEnd]

def m_point_2(v1,v2,p1,p2,t):
	"v1-v2 ist die Gerade, an der sich der Abstands-Punkt orientiert \
	der Abstands-Punkt liegt auf der Geraden p1,p2 \
	Rueckgabe ist: \
	stat = 1 wenn Punkt innerhalb von v1,v2 liegt \
	= 3 wenn Punkt parallel laeuft ausserhalb der Distanz \
	mm = Verhaltniss auf p1,p2 \
	p = der Punkt auf der Geraden"
	global me
	# Stat ist genauso wie stat bei mLine
	stat = 3
	mm = -10
	
	k0 = distAA2Edge(p1, v1, v2)
	k1 = distAA2Edge(p2, v1, v2)
	# m0, m1 Lotdistanz der beiden Punkte
	m0 = k0[0]
	m1 = k1[0]
	
	if k0[2] and k1[2]:
		# Test ob die beiden Punkte des Linienradius liegen
		if k0[0] < t and k1[0] < t:
			stat = 2
	
	p = [0,0,0]
	
	#print
	print "Lotdistanzen", m0, m1
	
	if ((m0 <= t and t <= m1) or (m1 <= t and t <= m0)) and m1 != m0:
		# einzige Vorraussetzung fuer weiteren Test ist, dass einer der beiden Punkte
		# innerhalb und einer ausserhalb der Lotdistanz liegt
		stat = 0
		m = (t - m0) / (m1 - m0)
		mm = m
		p = propEdgeV(p1, p2, m)
			
		k = distAA2Edge(p, v1, v2)
			
		if k[2]:
			stat = 1
			#lineDo(me, k[3], p)
				
	return [stat,mm,p]

def newVert(me, x, y, z, mode = 0):
	nonexists = True
	erg = 0 
	if mode == 0:
		for v in me.verts:
			if abs(x-v.co.x) < 0.001 and abs(y-v.co.y) < 0.001 and abs(z-v.co.z) < 0.001:
				erg = v
				nonexists = False
	if nonexists:
		erg = NMesh.Vert(x,y,z)
		me.verts.append(erg)	
	return erg
	
##########################
# 
##########################

# Objekte & Routinen

class mFace:
	def __init__(self,face):
		self.face = face
		self.rpoints = []
		self.adj = []
		self.norm = [0.0, 0.0, 0.0]
		# Berechnen der Normalen
		mi = [0.0, 0.0, 0.0]
		for v in face.v:
			mi = vecAdd(mi, pointToArray(v))
		mi = scalarMul(mi, 1.0 / len(face.v))
		di = [0.0, 0.0, 0.0]
		for i in range(len(face.v)):
			ip = (i+1) % len(face.v)
			d1 = vecSub(pointToArray(face.v[i]),mi)
			d2 = vecSub(pointToArray(face.v[ip]),mi)
			d = crossProd(d1,d2)
			di = vecAdd(di, d)
		self.norm = norm(di)

class mPoint:
	def __init__(self,co = [0.0, 0.0, 0.0]):
		self.p = None
		self.co = co
		self.x = False
		self.y = False
		self.z = False
		self.rand = False
		self.next = []
		self.edges = []
		
	def point(self, p):
		self.p = p
		self.co = [p.co.x, p.co.y, p.co.z]

	def setType(self, x, y, z):
		self.x = x
		self.y = y
		self.z = z
		self.co = [p.co.x, p.co.y, p.co.z]

	def getPoint(self, me):
		if (self.p == None):
			self.p = newVert(me, self.co[0], self.co[1], self.co[2])
		return self.p

	def getCo(self, me):
		return self.co
		
	def type(self, typeValue):
		if typeValue == 0:
			self.x = True
		elif typeValue == 1:
			self.y = True
		elif typeValue == 2:
			self.z = True

	def compare(self, other):
		"Gibt die Uebereinstimmung an \
		0 = keine Ubereinstimmung \
		1 = eine Koordinate \
		2 = zwei Koordinaten \
		bei 3 Koordinaten muesste es sich um denselben Punkt handeln"
		erg = 0
		ex = False
		ey = False
		ez = False
		v = 0.0001
		#if self.x and other.x and abs(self.co[0] - other.co[0]) < v:
		#	erg += 1
		#if self.y and other.y and abs(self.co[1] - other.co[1]) < v:
		#	erg += 1
		#if self.z and other.z and abs(self.co[2] - other.co[2]) < v:
		#	erg += 1

		if (self.x or other.x) and abs(self.co[0] - other.co[0]) < v:
			erg += 1
			ex = True
		if (self.y or other.y) and abs(self.co[1] - other.co[1]) < v:
			erg += 1
			ey = True
		if (self.z or other.z) and abs(self.co[2] - other.co[2]) < v:
			erg += 1
			ez = True
	 
		#if abs(self.co[0] - other.co[0]) < v:
		#	erg += 1
		#if abs(self.co[1] - other.co[1]) < v:
		#	erg += 1
		#if abs(self.co[2] - other.co[2]) < v:
		#	erg += 1
		return [erg, [ex, ey, ez]]


	def compareA(self, other):
		"Gibt die Uebereinstimmung an \
		0 = keine Ubereinstimmung \
		1 = eine Koordinate \
		2 = zwei Koordinaten \
		bei 3 Koordinaten muesste es sich um denselben Punkt handeln"
		erg = 0
		ex = FALSE
		ey = FALSE
		ez = FALSE
		v = 0.0001

		if (self.x and other.x) and abs(self.co[0] - other.co[0]) < v:
			erg += 1
			ex = TRUE
		if (self.y and other.y) and abs(self.co[1] - other.co[1]) < v:
			erg += 1
			ey = TRUE
		if (self.z and other.z) and abs(self.co[2] - other.co[2]) < v:
			erg += 1
			ez = TRUE
	 
		return [erg, [ex, ey, ez]]

	def mNextPoints(self,me,vp):
		"gibt eine Liste der verbundenen Punkte zurueck"
		erg = []
		for ee in me.edges:
			if ee.v1 == self.p: 
				erg.append(ee.v2)
			if ee.v2 == self.p:
				erg.append(ee.v1)
		return erg
	
def findMPoint(ppMenge, p):
	# sucht den Punkt p in der mPoint-Menge pp
	erg = None
	for pm in ppMenge:
		if pm.p == p:
			erg = pm
	return erg

def shortestWay(a):
	best = 10000
	erg = [0,[]]
	for i in range(len(a)):
		e = a[i]
		r = a[:i] + a[i+1:]
		neb = shortestWayP(e,r,best)
		if neb[0] < best:
			best = neb[0]
			erg = [best, [e] + neb[1]]
	return erg

def shortestWayP(pp,a,bestws, de = 0):
	erg = [0.0,[]]
	high = 100000
	for i in range(len(a)):
		e = a[i]
		d = dist(vecSub(pp.co, e.co)) 
		if d < bestws:
			r = a[:i] + a[i+1:]
			neb = shortestWayP(e, r, bestws-d, de + 1)
			if neb[0] + d < high:
				#bestws = neb[0] + d
				high = neb[0] + d
				erg = [high, [e] + neb[1]]
		else:
			erg[0] = d
	
	#print "erg: " , de , len(erg[1]), erg[0]
	return erg
			
class mLine:
	"Linie"
	def __init__(self,v1,v2):
		self.va = v1
		self.vb = v2
		self.dir = [0,0,0]
		self.polys = []

def addmLine(a,p1,p2):
	erg = findmLine(a,p1,p2)
	if 0 != erg:
		return erg  
	erg = mLine(p1,p2)
	a.append(erg)
	return erg

def findmLine(a,p1,p2):
	for li in a:
		if (li.va == p1 and li.vb == p2) or (li.va == p2 and li.vb == p1):
			return li
	return 0

class isoLinie:
	"Isolinie"
	def __init__(self, v1, v2, initDir = [False,False,False]):
		self.va = v1
		self.vb = v2
		self.p = []
		self.initDir = initDir

	def calc(self, v, are):
		p1 = pointToArray(self.va)
		p2 = pointToArray(self.vb)
		#print "*************"
		#printArray(p1)
		#printArray(p2)
		#print
		qu = 0.001
		for j in range(3):
			#print "**" , j
			mm = ceil(min (p1[j], p2[j]) / v) * v
			ma = floor(max (p1[j], p2[j]) / v) * v
			if -0.01 < ma - mm and are[j]:
				k = 0
				while (mm+k <= ma):
					m = mm+k
					m1 = abs(m-p1[j])
					m2 = abs(m-p2[j])
					msum = m1 + m2
					if qu < m1 or qu < m2:
						px = (p1[0] * m2 + p2[0] * m1) / msum
						py = (p1[1] * m2 + p2[1] * m1) / msum
						pz = (p1[2] * m2 + p2[2] * m1) / msum
						#print m, msum
						#printArray([px,py,pz])
						#print
						##p = newVert(me, px, py, pz)
						mp = mPoint([px,py,pz])
						#mp.point(p)
						mp.x = self.initDir[0]
						mp.y = self.initDir[1]
						mp.z = self.initDir[2]
						mp.type(j)
						self.p.append(mp)
					
					k += v
					
	def getAllIsoPoints(self):
		return self.p
					
	def drawLine(self, me):
		start = self.va
		startA = pointToArray(start)
		#aktuell = start
		erg = []
		b = True
		minH = 0
		while 0 < len(self.p):
			ib = 0
			best = 1000
			for i in range(len(self.p)):
				 hh = distAA(startA,self.p[i].co)
				 if hh < best:
					best = hh
					ib = i
			np = self.p.pop(ib)
			mayAdd = True
			if 0 < len(erg):
				if distAA(erg[-1].co,np.co) < 0.0001:
					mayAdd = False
			if mayAdd:
				erg.append(np)
		self.p = erg
		if 0 < len(self.p):
			next = self.p[0].getPoint(me)
			if self.va != next:
				me.addEdge(self.va, next)
			for i in range(len(self.p)-1):
				nnext = self.p[i+1].getPoint(me)
				if next != nnext:
					me.addEdge(next, nnext)
				next = nnext
			if self.vb != next:
				me.addEdge(next, self.vb)
		else:
			me.addEdge(self.va, self.vb)
			#aktuellA = pointToArray(aktuell)
			#b = False
			#h = 1000
			#for i in range(len(self.p)):
			#	hh = distAA(startA,self.p[i].co)
			#	if hh < h and minH + 0.0001 < hh:
			#		# es wurde ein weiterer Punkt gefunden
			#		h = hh
			#		ti = i
			#		b = True
			#if b:
			#	minH = h
			#	nextP = self.p[ti].getPoint(me)
			#	if start != nextP:
			#		me.addEdge(aktuell, nextP)
			#	aktuell = nextP
		#if aktuell != self.vb:
		#	me.addEdge(aktuell, self.vb)
			
	
def addLine(col, v1, v2):
	b = 0
	for m in col:
		if (m.va == v1 and m.vb == v2) or (m.va == v2 and m.vb == v1):
			b = 1
	if b == 0:
		mm = mLine(v1,v2)
		col.append(mm)
		
def addEdgeE(p0, p1, dist, mode = 0):
	# zieht eine Linie von p0 nach p1
	# hat dist einen Wert, so werden in festen Abstaenden dist
	# Zwischenpunkte gesetzt
	# mode hat (noch) keine Funktion
	global me
	if dist < 0.001:
		me.addEdge(p0, p1)
	else:
		a0 = pointToArray(p0)
		a1 = pointToArray(p1)
		da = distAA(a0,a1)
		seg = int(da / dist + 0.999)
		if mode == 0:
			v0 = p0
			if 0 < seg:
				m = scalarMul(norm(vecSub(a0,a1)), dist)
				for i in range(1, seg):
					av = vecAdd(a0, scalarMul(m,i))
					pp = NMesh.Vert(av[0], av[1], av[2])
					me.verts.append(pp)
					me.addEdge(v0, pp)
					v0 = pp 
			me.addEdge(v0, p1)
		else:
			v0 = p1
			if 0 < seg:
				m = scalarMul(norm(vecSub(a1,a0)), dist)
				for i in range(1, seg):
					av = vecAdd(a1, scalarMul(m,i))
					pp = NMesh.Vert(av[0], av[1], av[2])
					me.verts.append(pp)
					me.addEdge(v0, pp)
					v0 = pp 
			me.addEdge(v0, p0)
			
def isbetween(me,p1,pa,pb,p0):
	# testet, ob p1-p0 zwischen den beiden geraden p1-pa und p1-pb liegt
	rp1 = p1.getPoint(me)
	pass
			
def clockSort(ar):
	# Sortiert mPunkte im Uhrzeigersinn bzw Ggenuhrzeigersinn
	pass

def checkTorsion(f):
	# gibt den kleinsten Cosinus zwischen den Kreuzprodukt der Mittelpunktes
	# und zweier Ecken wieder...
	# wenn der Ruckgabewert 1 ist, dann ist das Face planar
	# Eine Wert nur nahe 1 ist bereits eine deutliche Torsion 
	vc = 1.0
	if len(f.v) == 4:	
		center = [0,0,0]
		dir = []
		for v in f.v:
			center = vecAdd(center, pointToArray(v))
		center = scalarMul(center, 1.0 / len(f.v))	  

		for i in range(len(f.v)):
			ip = (i+1) % len(f.v)
			a1 = pointToArray(f.v[i])
			a2 = pointToArray(f.v[ip])
			m = norm(crossProd(vecSub(center,a1),vecSub(center,a2)))
			dir.append(m)
	
		for i1 in range(len(dir)-1):
			for i2 in range(i1+1, len(dir)):
				v = scalarProd(dir[i1], dir[i2])
				if v < vc:
					vc = v
	return vc

def splitFace(me,f, mode = 0): 
	shortest = 10000
	shortestI = 0
	for i in range(len(f.v)):
		ip = (i+2) % len(f.v)
		a1 = pointToArray(f.v[i])
		a2 = pointToArray(f.v[ip])
		m = dist(vecSub(a2,a1))
		if m < shortest:
			shortest = m
			shortestI = i
			
	f1 = NMesh.Face()
	f2 = NMesh.Face()
	for i in range(3):
		i1 = (shortestI + i) %  len(f.v)
		i2 = (shortestI + i + 2) %  len(f.v)
		f1.v.append(f.v[i1])
		f2.v.append(f.v[i2])
	f1.materialIndex = f.materialIndex
	f2.materialIndex = f.materialIndex
	f1.mode = f.mode
	f2.mode = f.mode
	f1.smooth = f.smooth
	f2.smooth = f.smooth
	me.addFace(f1)
	me.addFace(f2)
	me.removeFace(f)
	return [f1,f2]
		  
###########################################
# Start Fill 2
###########################################

def isPartOfEdge(p1,p2,p):
	# ist Punkt p Teil der Kante (beschrieben durch p1, p2)
	d = [0, 0, 0]
	db = [False, False, False]
	qu = 0.001
	rqu = 0.00001
	for i in range(3):
		if qu < abs(p1[i] - p2[i]):
			d[i] = (p[i] - p2[i]) / (p1[i] - p2[i])
			db[i] = True
	stat = False
	# falls ein Wert nicht gesetzt werden konnte
	# dann nimmt er einen Wert einer anderen Stelle an
	for i in range(3):
		if not db[i]:
			for j in range(2):
				ip = (i+j+1) % 3
				# print i, j, ip
				if db[ip]:
					d[i] = d[ip]
					db[i] = True	
				
	inSquare = True
	for i in range(3):
		if p[i] < min(p1[i],p2[i]) - qu or max(p1[i],p2[i]) + qu < p[i]:
			inSquare = False
			
	if abs(d[0] - d[1]) < rqu and abs(d[0] - d[2]) < rqu and abs(d[1] - d[2]) < rqu \
	   and inSquare:
		stat = True
	#printArray(d)
	return stat  

def mpolyCenter(poly):
	m = [0,0,0]
	for p in poly:
		pa = pointToArray(p.p)
		m = vecAdd(m, pa)   
		
	m = scalarMul(m,1.0 / len(poly))
	return m

def drawMPoly(me, poly, mode = 0, tface = None):
	erg = []
	m = mpolyCenter(poly)
	
	if mode == 0:
		mm = newVert(me, m[0], m[1], m[2])
		for p in poly:
			me.addEdge(mm, p.p)
	elif mode == 1:
		mm = newVert(me, m[0], m[1], m[2])
		for pi in range(len(poly)):
			p1 = poly[pi]
			p2 = poly[(pi + 1) % len(poly)]
			pmm = scalarMul(vecAdd(p1.p,p2.p), 0.5)
			p = newVert(me, pmm[0], pmm[1], pmm[2])
			me.addEdge(mm, p)
	else:
		if len(poly) < 5:
			ff = NMesh.Face()
			for pi in range(len(poly)):
				ff.v.append(poly[pi].p)
			ff.smooth = 1
			if tface != None:
				ff.materialIndex = tface.materialIndex
				ff.mode = tface.mode
				ff.smooth = tface.smooth
			me.addFace(ff)
			erg.append(ff)
		else:
			dplus = 2
			if 5 < len(poly):
				dplus = 3
			boffs = 0
			dist = 0
			
			for pi in range(len(poly)):
				pe = (pi + dplus) % len(poly)
				d = distPP(poly[pi].p,poly[pe].p)
				if d > dist:
					boffs = pi
					dist = d
			boffsEnd = (dplus + boffs) % len(poly)
			#me.addEdge(poly[boffs].p, poly[(dplus + boffs) % len(poly)].p)
			
			ff = NMesh.Face()
			for pi in range(len(poly)):
				e = (pi + boffs) % len(poly)
				ff.v.append(poly[e].p)
				if (e == boffsEnd):
					fa = ff
					ff = NMesh.Face()
					ff.v.append(poly[e].p)
			ff.append(poly[boffs % len(poly)].p)
			fa.smooth = 1
			ff.smooth = 1
			if tface != None:
				fa.materialIndex = tface.materialIndex
				ff.materialIndex = tface.materialIndex
				fa.mode = tface.mode
				ff.mode = tface.mode
				fa.smooth = tface.smooth
				ff.smooth = tface.smooth
			me.addFace(fa)
			me.addFace(ff)
			erg.append(fa)
			erg.append(ff)
	return erg
				
def dirMLineToPoly(v1,v2,poly):
	"bestimmt die Richtung der Linie zum Polygon"
	global me
	m = mpolyCenter(poly)
	
	medPath = scalarMul(vecAdd(v1.p,v2.p),0.5)				  
	d1n = norm(vecSub(v2.p,v1.p))
	d2 = vecSub(medPath,m)
	vertLine = vecSub(scalarMul(d1n,scalarProd(d1n,d2)),m)
	d = vecSub(medPath,vertLine)
	#line (me,medPath,vecAdd(medPath,norm(d)))
	
	return d



def startDrawTest2():
	global me
	Blender.Window.EditMode(0)

	if me == 0:
		scn = Scene.GetCurrent()
		ob = scn.getActiveObject()
		me = ob.getData()
		Mesh.Mode(0)

	poi = []
	for p in me.verts:
		if p.sel == 1:
			poi.append(p)

	f = mFace(None)
	f.rpoints = poi
			
	gfill(me,f)
	
	me.update(1)
	Blender.Window.EditMode(1)
	me = 0


def gfill(me,poi):
	#global me
	
	#print "*** startFill ***"
	time_start = sys.time()
	
	# rpoi sind die auf Duplikate getesten Punkte
	rpoi = []
	for p in poi.rpoints:
		if p not in rpoi:
			rpoi.append(p)
	
	# Kante suchen
	# vp = Kanten
	# vpa = Kanten mLines
	vp = []
	vpa = []
	mPointListe = []
	anzahlKanten = 0
	
	#print len(poi.rpoints), len(rpoi)
	for p in rpoi:
		mPoint_p = findMPoint(mPointListe, p)
		if mPoint_p == None:
			mPoint_p = mPoint()
			mPoint_p.point(p)  
			mPointListe.append(mPoint_p)
		pa = getConnected(me, p)
		for pn in pa:
			if 0 < rpoi.count(pn):
				if p.index < pn.index:
					vp.append([p,pn])
					mayAdd = True
					## Kante darf nicht Teil des Originalfaces sein
					#for vi in range(len(poi.face.v)):
					#   vip = (vi + 1) % len(poi.face.v)
					#   if (poi.face.v[vi] == p and poi.face.v[vip] == pn) or (poi.face.v[vi] == pn and poi.face.v[vip] == p):
					#	   mayAdd = False
					
					# Es darf kein Punkt auf dieser Kante liegen
					for ptest in rpoi:
						if ptest != p and ptest != pn:
							#print pointPrint(p),
							#print pointPrint(pn),
							#print pointPrint(ptest),
							if isPartOfEdge(p,pn,ptest):
								mayAdd = False
								#pass
								#print mayAdd,
							#print
								
					mPoint_pn = findMPoint(mPointListe, pn)
					if mPoint_pn == None:
						mPoint_pn = mPoint()
						mPoint_pn.point(pn)  
						mPointListe.append(mPoint_pn)
						
					if mayAdd:
						anzahlKanten += 1
						#k = anzahlKanten / 10.0
						#line(me, [p[0]+k,p[1],p[2]], [pn[0]+k,pn[1],pn[2]])
						mLine = addmLine(vpa,mPoint_p,mPoint_pn)
						mPoint_p.next.append(mPoint_pn)
						mPoint_pn.next.append(mPoint_p)
						mPoint_p.edges.append(mLine)
						mPoint_pn.edges.append(mLine)
						
	# print anzahlKanten
	# print vpa
	
	if anzahlKanten < 2:
		vpa = [] 
		
	#for p in rpoi:
	#   print pointPrint(p)
	#   print p
						
	# von den mPoints die Nachbarmenge finden
	#for i1 in range(len(mPointListe)-1):
	#   for i2 in range(i1+1,len(mPointListe)):
	#	   p1 =  mPointListe[i1]
	#	   p2 =  mPointListe[i2]
	#	   for ee in me.edges:
	#		   if (ee.v1 == p1.p and ee.v2 == p2.p) or (ee.v2 == p1.p and ee.v1 == p2.p):
	#			   p1.next.append(p2)
	#			   p2.next.append(p1)
						
	#print "Struktur erfasst : %.2f"  % (sys.time() - time_start)
		
	allPolys = []
	
	# je hoher der Wert, umso spitzer duerfen die Kanten werten (aber irgendwas stimmt da noch nicht)
	edgeQual = 0.9
	
	#vpa = vpa[0:1]
	
	# Es wird nach und nach die Suchtiefe fur die
	# Polygone erweitert
	# jmm ist dabei die Anzahl der Kanten
	# jmb ist die Richtung
	for jm in range(0,8):
		jmm = 2 + jm / 2
		jmb = jm % 2
		# print jmm, jmb
		for vi in range(len(vpa)):
			v = vpa[vi]
			#print jmm, jmb, "******* Anzahl Polis" , len(v.polys) 
			if len(v.polys) < 2:
				for j in range(2 - len(v.polys)):
					pp = v.va
					pk = v.vb
					if jmb == 1:
						pp = v.vb
						pk = v.va
					# print pp,pk
					#pointPrint(pp.p)
					#pointPrint(pk.p)
					
					#print
					
					path = [pk,pp]
					touchedLines = [v]
							
					p = pp
					k = 0
					final = False
					while p != pk and k < jmm:
						k += 1
						# alle Nachbarn
						#paa = mp_getConnected(p)
						paa = p.next
						#print paa
						# pa sind die gefilterten Nachbarn
						pa = []
						for pj in paa:
							mayAdd = True
							for ppath in path:
								# Test. ob man den Pfad zuruckgeht
								if ppath == pj:
									mayAdd = False
									
							if mayAdd:
								# Test auf Manifold
								#vp = findmLine(vpa,p,pj)
								vp = findmLine(p.edges,p,pj)
								#sprint vp, len(vp.polys) 
								if 1 < len(vp.polys):
									mayAdd = False
									
							if mayAdd and k == 1 and 0.01 < dist(v.dir):
								# Test auf gleiche Richtung
								# die Qualitaet der Kante wird gemessen
								mpoly = path + [pj] 
								#print "!!!! Test Ausrichtung Flaeche"
								#printArray (v.dir)
								d = dirMLineToPoly(pp,pk,mpoly)
								#printArray (d)
								#print
								vc = vecCos(d,v.dir)
								if edgeQual < vc:
									mayAdd = False
								
							if mayAdd:
								pa.append(pj)
						
						# testen der besten Verbindung
						best = -2
						besthh = 10000
						pb = 0	  
						for pj in pa:
							#d = distPP(pp,pj) + distPP(pk,pj) 
							c = cosBetweenPoints(pp.p, pk.p, pj.p)
							if best < c:
								best = c
								pb = pj
							#if d < besthh:
							#   besthh = d
							#   pb = pj
							
						if pb != 0:
							#print "+ ", c, pointPrint(pp), pointPrint(pk), pointPrint(pj)
							
							# Sinus, um die Richtung zuzuordnen
							#d = sinBetweenPoints(path[-1], path[-2], pb)				   
							#lineDir(me, path[-1], norm(d))
												
							#print "=======", j, ", Pfad:", len(path), ", Knoten :" , pb.index,
							#printArray(v.dir)
							#print  

							p = pb
							path.append(p)
							if path[0] == pk:
								path = path[1:]
						else:
							p = 0
							k = 1000
						
				
					if p == pk:
						#print "Ziel erfolgreich erreicht"
						final = True 
					else:
						#print "Verlaufen"
						pass
				
					if final:
						# Pfad war erfolgreich, letzter Test, ob Poly sich deckt mit bisherigem Feld
						mayAdd = True
						#print "dirLineToPoly",pp,pk,path
						d = dirMLineToPoly(pp,pk,path)
						#printArray (d)
						#print
						vc = vecCos(d,v.dir)
						if edgeQual < vc:
							mayAdd = False
						
						if mayAdd:
							
							# Poly war erfolgreich
							allPolys.append(path)
							#vp = findmLine(vpa,v.va,v.vb)
							#vp.polys.append(path)
							#print medPath , d1n, d2
					
							# fur alle beteiligten Kanten das Poly anhangen
							for i in range(len(path)):
								ip = (i + 1) % len(path)
								v1 = path[i]
								v2 = path[ip]
						
								# Richtung der Flache bestimmen,
								# aufgrund der Richtung sind die weiteren
								# Flachen am Rand eingeschrankt
								d = dirMLineToPoly(v1,v2,path)
							
								#vp = findmLine(vpa,v1,v2)
								vp = findmLine(v1.edges,v1,v2)
								if 0 != vp: 
									vp.dir = d
									vp.polys.append(path)
									
		#print "Loop:",jmm,"Kanten, %.2f"  % (sys.time() - time_start)
		
	erg = []
	for i in range(len(allPolys)):
		poly = allPolys[i]
		erg.extend(drawMPoly(me, poly, 2, poi.face))
	return erg
			
	#print "finished : %.2f"  % (sys.time() - time_start)
		
###########################################
# Start Routine
###########################################

def iso(me,selfaces,stv,dirV,mode = 0):  
	isoV = [dirV[0] == 1, dirV[1] == 1, dirV[2] == 1]
	faces = []
	for f in selfaces:
		ct = checkTorsion(f)
		if ct < 0.99999 and mode == 0:
			ff = splitFace(me,f)
			faces.append(mFace(ff[0]))
			faces.append(mFace(ff[1]))
		else:
			faces.append(mFace(f))
			
	for f in faces:
		isolinien = []
		#print "***** Step 1 *****"
		for i in range(len(f.face.v)):
			ip = (i+1) % len(f.face.v)
			p1 = pointToArray(f.face.v[i])
			p2 = pointToArray(f.face.v[ip])
			
			il = isoLinie(f.face.v[i], f.face.v[ip], [False,False,False])
			il.calc(stv , isoV)
			isolinien.append(il)
			
		#print "***** Step 2 *****"
		
		# fp sind die Isopunkte der ersten Ordnung
		# fpo sind die Isopunkte der zweiten Ordnung
		#	2.te Ordnung = 2 Koordinaten stimmen ueberein
		fp = []
		fpo = []
		iso2 = []
		for iso in isolinien:
			fp.extend(iso.getAllIsoPoints())
		for i1 in range(len(fp)-1):
			for i2 in range(i1+1, len(fp)):
				c = fp[i1].compare(fp[i2])
				if (0 < c[0]):
					p1 = fp[i1].getPoint(me)
					p2 = fp[i2].getPoint(me)
					if p1 != p2:
						
						#me.addEdge(p1, p2)
					
						il = isoLinie(p1, p2,c[1])
						il.calc(stv , isoV)
						iso2.append(il)
		
		for iso in iso2:
			fpo.extend(iso.getAllIsoPoints())

		# Berechnen der gemittelten Punkte
		#for i1 in range(len(fpo)-1):
		#	for i2 in range(i1+1, len(fpo)):
		#		c = fpo[i1].compareA(fpo[i2])
		#		#print "***"
		#		#printArray(fpo[i1].co)
		#		#print [fpo[i1].x, fpo[i1].y, fpo[i1].z] , ",",
		#		#printArray(fpo[i2].co)
		#		#print [fpo[i2].x, fpo[i2].y, fpo[i2].z] , " = ",
		#		#print c[0]
		#		
		#		#if (1 < c[0]):
		#		#   p1 = fpo[i1].getCo(me)
		#		#   p2 = fpo[i2].getCo(me)
		#		#   pm = scalarMul(vecAdd(p1,p2),0.5)
		#		#   #printArray(pm)
		#		#   #print  
		#		#   ppm = newVert(me, pm[0], pm[1], pm[2])
		#		#   fpo[i1].point(ppm)
		#		#   fpo[i2].point(ppm)
		
		#print "***** Step 3 *****"
		
		# Die Linien der Kanten zeichnen
		for iso in isolinien:
			iso.drawLine(me)
					
		# Die Linien auf der Flaeche zeichnen
		for iso in iso2:
			iso.drawLine(me)
	
		# Alle Punkte sammeln
		#punkte = []
		for fppu in f.face.v:
			f.rpoints.append(fppu)
		for fpp in fp:
			fpp.rand = True
			f.rpoints.append(fpp.getPoint(me))
		for fpp in fpo:
			fpp.rand = False
			f.rpoints.append(fpp.getPoint(me))
			#if fpp.p != None:
			#   punkte.append(fpp.p)
			
		#ergPoints.append(punkte)
		#gfill(me,punkte)
		#print len(punkte)

	#print "*** Ausfuellen der Gitter ***" 
	#print "time now : %.2f"  % (sys.time() - time_start) 
	me.update(1)

	allNewFaces = []
		 
	for f in faces:
		allNewFaces.extend(gfill(me,f))

	#print "***** Step 4 loeschen der Faces *****"		  
	for f in faces:
		me.removeFace(f.face)

	return allNewFaces
		
me = 0
	
def startDraw(stv,dirV,mode):
	global me
	print "*** start ***"
	time_start = sys.time()
	Blender.Window.EditMode(0)
	
	if me == 0:
		scn = Scene.GetCurrent()
		ob = scn.objects.active
		me = ob.getData()
	
	Mesh.Mode(1)
		
	ergPoints = []

	selFaces = []
	for f in me.faces:
		if f.sel:
			selFaces.append(f)

	if mode == 1:
		selFaces = iso(me,selFaces,stv,dirV,0)
	elif mode == 2 or mode == 3:
		f1 = []
		f2 = []
		for f in selFaces:
			ct = checkTorsion(f)
			if ct < 0.99999:
				f1.append(f)
			else:
				f2.append(f)
		f2 = iso(me,f2,stv,dirV,1)
		if mode == 2:
			f1 = iso(me,f1,stv,[0,0,dirV[2]],1)
			f1 = iso(me,f1,stv,[0,dirV[1],0],1)
			f1 = iso(me,f1,stv,[dirV[0],0,0],1)
		if mode == 3:
			f1 = iso(me,f1,stv,[0,dirV[1],0],1)
			f1 = iso(me,f1,stv,[dirV[0],0,0],1)
			f1 = iso(me,f1,stv,[0,0,dirV[2]],1)
		selFaces = f1 + f2

	me.update(1)
	for f in selFaces:
		for v in f.v:
			v.sel = 1
		
	me.update(1)
	Blender.Window.EditMode(1)
	me = 0
	print "finished : %.2f"  % (sys.time() - time_start)
	
def markFaces():	
	global me
	print "*** start ***"
	time_start = sys.time()
	Blender.Window.EditMode(0)
	
	if me == 0:
		scn = Scene.GetCurrent()
		ob = scn.objects.active
		me = ob.getData()
	
	Mesh.Mode(1)

	allFaces = []
	for f in me.faces:
		allFaces.append(mFace(f))
		
	# eine Zwischenspeicherung aller Punkte und der angrenzenden mFaces
	allPoints = []
	for p in me.verts:
		pp = [p,[]]
		allPoints.append(pp)
		
	for f in allFaces:
		ff = f.face
		for v in ff.v:
			for ii in range(len(allPoints)):
				if allPoints[ii][0] == v:
					allPoints[ii][1].append(f)
					
	#for v in allPoints:
	#	print "*****"
	#	for f in v[1]:
	#		print f
					
	# Nachbarn setzen
	for f in allFaces:
		#print "****"
		ff = f.face
		for i in range(len(ff.v)):
			ip = (i+1) % len(ff.v)
			# v1 und v2 sind die Punkte einer Kante
			# vp1 und vp2 sind die Faces, die an einen Punkt angrenzen
			v1 = ff.v[i]
			v2 = ff.v[ip]
			vp1 = None
			vp2 = None
			for vv in allPoints:
				if vv[0] == v1:
					vp1 = vv[1]
				if vv[0] == v2:
					vp2 = vv[1]
			if vp1 != None and vp2 != None:
				nFace = 0
				for f1 in vp1:
					for f2 in vp2:
						if f1 == f2 and f1 != f:
							nFace = f1
				if nFace != 0:
					#print nFace
					f.adj.append(nFace)
		#print (len(f.adj))
					
	#
	b = True
	while b:
		b = False
		for f in allFaces:
			#print f, len(f.adj)
			if f.face.sel == 1:
				ff = f.face
				for fn in f.adj:
					if fn.face.sel == 0:
						r = abs(scalarProd(f.norm,fn.norm))
						if 0.9999 < r:
							fn.face.sel = 1
							b = True
							for v in fn.face.v:
								v.sel = 1
			
				

	Mesh.Mode(4)
	me.update(1)
	Blender.Window.EditMode(1)
	me = 0
	print "finished : %.2f"  % (sys.time() - time_start)
		
#******************************

class fLine:
	"Linie"
	def __init__(self,v1,v2):
		self.va = v1
		self.vb = v2
		self.faces = []
		self.adjEdges = []
		self.marked = False
		
	def isEdge(self):
		m = 0
		for f in self.faces:
			if f.sel == 1:
				m += 1
		return m < 2

def addfLine(a,p1,p2):
	erg = findfLine(a,p1,p2)
	if 0 != erg:
		return erg  
	erg = fLine(p1,p2)
	a.append(erg)
	return erg

def findfLine(a,p1,p2 = 0):
	if p2 == 0:
		# alle Linien, die den Punkt p1 enthalten
		erg = []
		for li in a:
			if li.va == p1 or li.vb == p1:
				erg.append(li)
		return erg
	else:
		for li in a:
			if (li.va == p1 and li.vb == p2) or (li.va == p2 and li.vb == p1):
				return li
	return 0

def reduceFaces():	
	global me
	print "*** start ***"
	time_start = sys.time()
	Blender.Window.EditMode(0)
	
	if me == 0:
		scn = Scene.GetCurrent()
		ob = scn.objects.active
		me = ob.getData()
	
	Mesh.Mode(1)
	
	allEdges = []
	allMarkedFaces = []
	
	# Erstellen der Edges
	for f in me.faces:
		if f.sel == 1:
			allMarkedFaces.append(f)
			for i in range(len(f.v)):
				ip = (i+1) % len(f.v)
				nl = addfLine(allEdges, f.v[i], f.v[ip])
				nl.faces.append(f)
	
				
	# Erstellen der Kantenzuege
	allRounds = []
	for l in allEdges:
		#print l
		if not l.marked:
			if l.isEdge():
				newRound = [l]
				stopPoint = l.va
				next = l.vb
				lb = l
				c = 0
				while next != stopPoint and c<100:
					c += 1
					# alle angrenzenden Kanten
					ll = findfLine(allEdges, next)
					ln = 0
					for lnn in ll:
						# die naechste Kante ist auch am Rand und
						# darf nicht die letzte Kante sein
						if lnn != lb and lnn.isEdge():
							ln = lnn 
					if ln == 0:
						next = stopPoint
					else:
						lb = ln
						if next == lb.va:
							next = lb.vb
						else:
							next = lb.va
						newRound.append(lb)
						lb.marked = True
				allRounds.append(newRound)
				
	#print allRounds
	
	# bei den Kantenzuegen die Eckpunkte feststellen
	edgeRounds = []
	for r in allRounds:
		edgeR = []
		for i in range(len(r)):
			ip = (i+1) % len(r)
			l1 = r[i]
			l2 = r[ip]
			d1 = vecSub(pointToArray(l1.va), pointToArray(l1.vb))
			d2 = vecSub(pointToArray(l2.va), pointToArray(l2.vb))
			c = abs(vecCos(d1,d2))
			if c < 0.99999:
				p = 0
				if l1.va == l2.va:
					p = l1.va
				elif l1.vb == l2.vb:
					p = l1.vb
				elif l1.vb == l2.va:
					p = l1.vb
				else:
					p = l1.va
				if p != 0:
					edgeR.append(pointToArray(p))
		edgeRounds.append(edgeR)
	
	Mesh.Mode(1)
	
	# alle Faces loeschen
	# komischerweise werden dabei nicht die freiwerdenden Punkte geloescht
	allDPoints = []
	for f in allMarkedFaces:
		for v in f.v:
			if allDPoints.count(v) == 0:
				allDPoints.append(v)
		me.removeFace(f)
	
	# alle Verts aus der Liste entfernen, die noch irgendwie an einem Face haengen
	for f in me.faces:
		for v in f.v:
			if 0 < allDPoints.count(v):
				allDPoints.remove(v)
				
	for v in allDPoints:
		me.verts.remove(v)
	
		
	for v in me.verts:
		v.sel = 0 
					
	isFillable = False
	
	if len(edgeRounds) == 1:
		if len(edgeRounds[0]) < 5:
			isFillable = True
					
	for r in edgeRounds:
		pp = []
		for co in r:
			np = NMesh.Vert(co[0], co[1], co[2])
			np.sel = 1
			me.verts.append(np)	
			pp.append(np)
		if isFillable:		
			ff = NMesh.Face()
			for p in pp:
				ff.v.append(p)
			ff.smooth = 1
			me.addFace(ff)
		else:
			for i in range(len(pp)):
				ip = (i+1) % len(pp)			
				me.addEdge(pp[i], pp[ip])
				
	
	#Mesh.Mode(4)
			
	me.update(1)
	Blender.Window.EditMode(1)
	me = 0
	print "finished : %.2f"  % (sys.time() - time_start)

#Draw.Register 
Register(draw, event, bevent)