[Bf-committers] Surface and Volume modeling with new developed Bezier-Surface math

Sun Sep 29 18:46:14 CEST 2013

```  Hello,

I am nowhere with the web page that I wanted to create. However I
started to dig a bit into blender scripting and so I created a script
which makes the thing more visible.

The script creates tangential planes and handles on the vertices of
the selected mesh.  However it does not create bezier surfaces or volume
yet.

The name is: Bezier Volume Operator

Sry if the code may look funny. I never programmed phyton before nor I
know what blender offers (math functions and data structures).

regards

--------------------

import bpy
import math

class XYZPoint:

def __init__(self,x,y,z):
self.x = x
self.y = y
self.z = z

description = "XYZPoint stores the x,y and z coord of a point or a
vector"

class HNFPlane:

def __init__(self,a,b,c,d):
self.a = a
self.b = b
self.c = c
self.d = d

description = "HNFPlane a*x + b*y + c* z + d = 0, whereas
a*a+b*b+c*c =1"

def calcdistance(self,x,y,z):
dist = self.a * x + self.b * y + self.c * z + self.d
return dist

def flipnormal(self):
self.a = -self.a
self.b = -self.b
self.c = -self.c
self.d = -self.d

def flipnormaltowardpoint(self, x,y,z):
dist = self.calcdistance(x,y,z)
if dist < 0:
self.flipnormal()

def moveplanetopoint(self, x,y,z):
dist = self.calcdistance(x,y,z)
newd = self.d - dist * (self.a*self.a + self.b*self.b +
self.c*self.c)
self.d = newd

def CalculateXYZGlobalAngles(self):
vBackTurnAngleY = - math.atan2(self.a, self.c)
vNewZ   =   math.sqrt(self.a * self.a + self.c * self.c);
vBackTurnAngleX =   math.atan2(self.b, vNewZ)
GlobalAngleX = - vBackTurnAngleX;
GlobalAngleY = - vBackTurnAngleY;
return (GlobalAngleX, GlobalAngleY, 0)

class BezierVolumeOperator(bpy.types.Operator):
bl_idname = "object.beziervolume_operator"
bl_label = "Bezier Volume Operator"

def execute(self, context):
print("Start")
HandleLengthFactor = 0.333 #Factor of edge lengh
SourceData = context.active_object.data
VertIndexOfInterest = 0;
for VertOfInterest in SourceData.vertices:
TangentialPlaneVerts  =
ExtractTangentialVerticesArray(SourceData, VertIndexOfInterest)
TangentialHNFPlane    = CreateTangentialPlane(VertOfInterest,
TangentialPlaneVerts)
HandelsVectors        =
CalculateHandleVectors(TangentialHNFPlane, VertOfInterest,
TangentialPlaneVerts, HandleLengthFactor)
drawplane(context, VertOfInterest, TangentialHNFPlane,
HandelsVectors)
#print(TangentialPlaneVerts)
VertIndexOfInterest += 1
print("Stop")
return {'FINISHED'}

def drawplane(context, destver, iHNFPlane, iHandelsVectors):
vGlobalPlaneAngles = iHNFPlane.CalculateXYZGlobalAngles();

#rotate the handles to global plane
globalizedHhandels = []
for handle in iHandelsVectors:
globalhandle = TurnToGlobal(handle, vGlobalPlaneAngles)
globalizedHhandels.append(globalhandle)

maxhandlelength = 0
for gbh in globalizedHhandels :
length = math.sqrt(gbh[0] * gbh[0] + gbh[1] * gbh[1] +
gbh[2] * gbh[2])
if  length > maxhandlelength:
maxhandlelength = length

coords = []
vAngles = (0,45,90,135,180,225,270,315)
for vAngle in vAngles:
vX = radius * math.cos(3.1415926 * vAngle/180)
vY = radius * math.sin(3.1415926 * vAngle/180)
coords.append((vX,vY,0))
#coords = [(-1,-1,0), (1,-1,0), (1,1,0), (-1,1,0)]
faces = [(0,1,2,3,4,5,6,7)]

vertexcount = 8
vtrianglebasefactor = 0.1
for globhandle in globalizedHhandels:
point1 = ( vtrianglebasefactor*globhandle[1], -
vtrianglebasefactor*globhandle[0], 0)
point2 = (-vtrianglebasefactor*globhandle[1],
vtrianglebasefactor*globhandle[0], 0)
coords.append(point1)
coords.append(globhandle)
coords.append(point2)
faces.append((vertexcount,vertexcount+1,vertexcount+2))
vertexcount += 3
#print(coords)
#print(faces)
me = bpy.data.meshes.new("TangPlaneMesh")
ob = bpy.data.objects.new("TangPlaneObject", me)
ob.rotation_euler = vGlobalPlaneAngles
ob.location = (destver.co[0],destver.co[1],destver.co[2])
ob.parent = context.active_object
me.from_pydata(coords,[],faces)
me.update(calc_edges=True)
return {'FINISHED'}

def TurnToGlobal(handle, vGlobalPlaneAngles):
vX = handle[0]
vY = handle[1]
vZ = handle[2]
# Z first
vAlfaSin = math.sin(-vGlobalPlaneAngles[2]);
vAlfaCos = math.cos(-vGlobalPlaneAngles[2]);
vNewX = vAlfaCos * vX - vAlfaSin * vY;
vNewY = vAlfaSin * vX + vAlfaCos * vY;
vX = vNewX
vY = vNewY
# Y first
vAlfaSin = math.sin(-vGlobalPlaneAngles[1]);
vAlfaCos = math.cos(-vGlobalPlaneAngles[1]);
vNewZ = vAlfaCos * vZ - vAlfaSin * vX;
vNewX = vAlfaSin * vZ + vAlfaCos * vX;
vZ = vNewZ
vX = vNewX
# X after
vAlfaSin = math.sin(-vGlobalPlaneAngles[0]);
vAlfaCos = math.cos(-vGlobalPlaneAngles[0]);
vNewY = vAlfaCos * vY - vAlfaSin * vZ;
vNewZ = vAlfaSin * vY + vAlfaCos * vZ;
vY = vNewY
vZ = vNewZ
return (vX,vY,vZ)

def ExtractTangentialVerticesArray(SourceData, indexofinterest):
verticesindex = []
arrvertices = []
for ed in SourceData.edges:
found = False
for vert in ed.vertices:
if   indexofinterest == vert:
found = True
if (found):
for vert in ed.vertices:
if   indexofinterest != vert:
verticesindex.append(vert)
for vert in verticesindex:
arrvertices.append(SourceData.vertices[vert])
return arrvertices

def CreateTangentialPlane(destver, planevertarray):
cEpsilon = 0.00000001

arraysize = 0.0
for vert in planevertarray:
arraysize += 1.0

if arraysize < cEpsilon:
return HNFPlane(0,0,1,0)

vDN = arraysize
vSx = 0.0
vSy = 0.0
vSx2 = 0.0
vSy2 = 0.0
vSz2 = 0.0
vSxy = 0.0
vSz = 0.0
vSxz = 0.0
vSyz = 0.0

for vert in planevertarray:
x = vert.co[0]
y = vert.co[1]
z = vert.co[2]
vSx  += x
vSy  += y
vSz  += z
vSx2 += x*x
vSy2 += y*y
vSz2 += z*z
vSxy += x*y
vSxz += x*z
vSyz += y*z

#[a,b,c]
vNennerABC = (vSx2 * (vSyz * vSyz - vSy2 * vSz2) + vSxy * vSxy *
vSz2 - 2.0 * vSxy * vSxz * vSyz + vSxz * vSxz * vSy2)
#[b,c,d]
vNennerBCD = (vSy2 * (vSz * vSz - vSz2 * vDN) + vSyz * vSyz * vDN +
vSy * vSy * vSz2 - 2.0 * vSy * vSyz * vSz)
#[a,c,d]
vNennerACD = (vSx2 * (vSz * vSz - vSz2 * vDN) + vSxz * vSxz * vDN +
vSx * vSx * vSz2 - 2.0 * vSx * vSxz * vSz)
#[a,b,d]
vNennerABD = (vSx2 * (vSy * vSy - vSy2 * vDN) + vSxy * vSxy * vDN +
vSx * vSx * vSy2 - 2.0 * vSx * vSxy * vSy)

vMaxNenner = abs(vNennerABC)
if abs(vNennerBCD) > vMaxNenner:
vMaxNenner = abs(vNennerBCD)
if abs(vNennerACD) > vMaxNenner:
vMaxNenner = abs(vNennerACD)
if abs(vNennerABD) > vMaxNenner:
vMaxNenner = abs(vNennerABD)

if vMaxNenner< cEpsilon:
return HNFPlane(0,0,1,0)

vA = 0.0
vB = 0.0
vC = 0.0
vD = 0.0

if abs(vNennerABC) >= vMaxNenner:
vFactor = 1.0 / vNennerABC
vD = 1.0
vA = vFactor *  (vD * vSx  * (vSy2 * vSz2 - vSyz * vSyz) + vD *
vSxy * (vSyz * vSz - vSy * vSz2) + vD * vSxz * (vSy * vSyz - vSy2 * vSz))
vB = vFactor * -(vD * vSx2 * (vSyz * vSz - vSy * vSz2) + vD *
vSx * vSxy * vSz2 + vSxz * (-vD * vSxy * vSz - vD * vSx * vSyz) + vD *
vSxz * vSxz * vSy)
vC = vFactor *  (vD * vSx2 * (vSy2 * vSz - vSy * vSyz)   - vD *
vSxy * vSxy * vSz + vD * vSx * vSxy * vSyz + vSxz * (vD * vSxy * vSy -
vD * vSx * vSy2))

elif abs(vNennerBCD) >= vMaxNenner:
vFactor = 1.0 / vNennerBCD
vA = 1.0
vB = vFactor *  (vA * vSxy * (vSz2 * vDN - vSz * vSz) + vA * vSyz
* (vSx * vSz - vSxz * vDN) + vA * vSy * (vSxz * vSz - vSx * vSz2))
vC = vFactor *  (vA * vSy2 * (vSxz * vDN - vSx * vSz) + vSyz *
(vA * vSx * vSy - vA * vSxy * vDN) + vA * vSxy * vSy * vSz - vA * vSxz *
vSy * vSy)
vD = vFactor * -(vA * vSy2 * (vSxz * vSz - vSx * vSz2) + vA *
vSxy * vSy * vSz2 + vSyz * (-vA * vSxy * vSz - vA * vSxz * vSy) + vA *
vSx * vSyz * vSyz)

elif abs(vNennerACD) >= vMaxNenner:
vFactor = 1.0 / vNennerACD
vB = 1.0
vA = vFactor *  (vB * vSxy * (vSz2 * vDN - vSz * vSz) + vB * vSxz
* (vSy * vSz - vSyz * vDN) + vB * vSx * (vSyz * vSz - vSy * vSz2))
vC = vFactor * -(vB * vSx2 * (vSy * vSz - vSyz * vDN) + vSxz *
(vB * vSxy * vDN - vB * vSx * vSy) - vB * vSx * vSxy * vSz + vB * vSx *
vSx * vSyz)
vD = vFactor * -(vB * vSx2 * (vSyz * vSz - vSy * vSz2) + vB * vSx
* vSxy * vSz2 + vSxz * (-vB * vSxy * vSz - vB * vSx * vSyz) + vB * vSxz
* vSxz * vSy)

elif abs(vNennerABD) >= vMaxNenner:
vFactor = 1.0 / vNennerABD
vC = 1.0;
vA = vFactor * (vC * vSxy * (vSy * vSz - vSyz * vDN) + vC * vSxz
* (vSy2 * vDN - vSy * vSy) + vC * vSx * (vSy * vSyz - vSy2 * vSz))
vB = vFactor * -(vC * vSx2 * (vSy * vSz - vSyz * vDN) + vSxz *
(vC * vSxy * vDN - vC * vSx * vSy) - vC * vSx * vSxy * vSz + vC * vSx *
vSx * vSyz)
vD = vFactor * (vC * vSx2 * (vSy2 * vSz - vSy * vSyz) - vC * vSxy
* vSxy * vSz + vC * vSx * vSxy * vSyz + vSxz * (vC * vSxy * vSy - vC *
vSx * vSy2))

vNormalize = math.sqrt(vA * vA + vB * vB + vC * vC);
if vNormalize < cEpsilon:
return HNFPlane(0,0,1,0)

vNormalizeFactor = 1.0 / vNormalize
vA = vNormalizeFactor * vA
vB = vNormalizeFactor * vB
vC = vNormalizeFactor * vC
vD = vNormalizeFactor * vD

vHNFPlane = HNFPlane(vA,vB,vC,vD)

vHNFPlane.flipnormaltowardpoint(destver.co[0], destver.co[1],
destver.co[2])

vHNFPlane.moveplanetopoint(destver.co[0], destver.co[1],
destver.co[2])

#print(vHNFPlane.a)
#print(vHNFPlane.b)
#print(vHNFPlane.c)
#print(vHNFPlane.d)

return  vHNFPlane

def CalculateHandleVectors(iTangentialHNFPlane, iVertOfInterest,
iTangentialPlaneVerts, iHandleLengthFactor):
Handles = []
for vert in iTangentialPlaneVerts:
dist = iTangentialHNFPlane.calcdistance(vert.co[0], vert.co[1],
vert.co[2])
planehitpoint = (vert.co[0] - dist*iTangentialHNFPlane.a,
vert.co[1] - dist*iTangentialHNFPlane.b, vert.co[2] -
dist*iTangentialHNFPlane.c)
handlevector  = (planehitpoint[0] - iVertOfInterest.co[0],
planehitpoint[1] - iVertOfInterest.co[1], planehitpoint[2] -
iVertOfInterest.co[2])
reducedhandle = (iHandleLengthFactor*handlevector[0],
iHandleLengthFactor*handlevector[1], iHandleLengthFactor*handlevector[2])
Handles.append(reducedhandle)
return Handles

def register():
bpy.utils.register_class(BezierVolumeOperator)

if __name__ == "__main__":
register()

```