nortikin / sverchok

Sverchok
http://nortikin.github.io/sverchok/
GNU General Public License v3.0
2.21k stars 234 forks source link

Alternative Recalc Normals sudocode #904

Closed zeffii closed 7 years ago

zeffii commented 7 years ago

should be a mesh util,

properties.

algorithm

  for each mesh:
    get first normal found
    compare each following normal, reverse order of verts found if normal does not broadly compare
  for each mesh
    get normal of polygon, if the average vector generated by connecting each vertex of the 
    polygon to the nearest attractor matches the found normal then keep index order, else reverse

   ( or some variation of this )

here the red dot is a single attractor

image

multiple attractors per mesh

image

zeffii commented 7 years ago

some manual normal finding, instead of using the Vector class..

here restated a few times to allow some speed comparissons.

def obtain_normal(p1, p2, p3):
    # http://stackoverflow.com/a/8135330/1243487
    x1, y1, z1 = p1
    x2, y2, z2 = p2
    x3, y3, z3 = p3
    nx = (y2 - y1)*(z3 - z1) - (z2 - z1)*(y3 - y1)
    ny = (z2 - z1)*(x3 - x1) - (x2 - x1)*(z3 - z1)
    nz = (x2 - x1)*(y3 - y1) - (y2 - y1)*(x3 - x1)
    return [nx, ny, nz]

def obtain_normal2(p1, p2, p3):
    # http://stackoverflow.com/a/8135330/1243487
    x2x1 = p2[0] - p1[0]
    x3x1 = p3[0] - p1[0]
    y2y1 = p2[1] - p1[1]
    y3y1 = p3[1] - p1[1]
    z2z1 = p2[2] - p1[2]
    z3z1 = p3[2] - p1[2]
    nx = y2y1*z3z1 - z2z1*y3y1
    ny = z2z1*x3x1 - x2x1*z3z1
    nz = x2x1*y3y1 - y2y1*x3x1
    return [nx, ny, nz]

def obtain_normal3(p1, p2, p3):
    # http://stackoverflow.com/a/8135330/1243487
    return [
        ((p2[1]-p1[1])*(p3[2]-p1[2]))-((p2[2]-p1[2])*(p3[1]-p1[1])),
        ((p2[2]-p1[2])*(p3[0]-p1[0]))-((p2[0]-p1[0])*(p3[2]-p1[2])),
        ((p2[0]-p1[0])*(p3[1]-p1[1]))-((p2[1]-p1[1])*(p3[0]-p1[0]))
    ]
nortikin commented 7 years ago

i see great idea with attractor points!

zeffii commented 7 years ago

small bit of test code to confirm normals calc is OK:
https://gist.github.com/7509dbe6293ad4ebb996f02b7237c1d0

the obtain_normal* code doesn't auto normalize, which might be interesting as a feature, the normal direction are correct, but the length is a function of the polygon surface area

def normalize(v):
    l = math.sqrt((v[0] * v[0]) + (v[1] * v[1]) + (v[2] * v[2]))
    return [v[0]/l, v[1]/l, v[2]/l]
zeffii commented 7 years ago

just sudocode for now

def perform_recalcs(input_data):
    num_data_channels = len(input_data)
    func = {3: homogenous, 4: attracted}.get(num_data_channels)
    return func(*input_data) if func else []

def homogenous(vlist, plist, flipped):
    '''
    - recalculate homogenous

    assume all normals are perpendicular to a single plane, but some normals 
    are flipped. This will turn them all to face the direction of the first found
    normal (or reverse if flipped is True).

    useful for islands of potentially flipped normals, where direction can not be 
    usefully inferred by vacinity.
    '''
    for each mesh:
        get first normal found
        compare each following normal, reverse order of verts found if normal does not broadly compare
    return corrected_polygons

def attracted(vlist, plist, flipped, attractors):
    '''
    - recalculate polygon direction based on attractors

    The median of the polygon is used to establish the location of a polygon.
    Each attractor is added to a kdtree, and for each median it's closest attractor
    is found and this relationship is used to determine if the polygon's direction
    should be flipped or not.

    '''

    for each mesh
        n = get normal of polygon, 
        if the average vector generated by connecting each vertex of the 
        polygon to the nearest attractor matches the found normal then keep index order, 
        else reverse
zeffii commented 7 years ago

some movement..

def obtain_normal3(p1, p2, p3):
    return [
        ((p2[1]-p1[1])*(p3[2]-p1[2]))-((p2[2]-p1[2])*(p3[1]-p1[1])),
        ((p2[2]-p1[2])*(p3[0]-p1[0]))-((p2[0]-p1[0])*(p3[2]-p1[2])),
        ((p2[0]-p1[0])*(p3[1]-p1[1]))-((p2[1]-p1[1])*(p3[0]-p1[0]))
    ]

def perform_recalcs(input_data):
    num_data_channels = len(input_data)
    func = {3: homogenous, 4: attracted}.get(num_data_channels)
    return func(*input_data) if func else []

def is_reasonably_opposite(n, normal_one):
    ...

def homogenous(vlist, plist, flipped):
    '''
    - recalculate homogenous

    assume all normals are perpendicular to a single plane, but some normals 
    are flipped. This will turn them all to face the direction of the first found
    normal (or reverse if flipped is True).

    useful for islands of potentially flipped normals, where direction can not be 
    usefully inferred by vacinity.

    unfortunately, this doesn't work great on concave polygons..
    '''

    h_polygons = []
    addp = h_polygons.append

    addp(plist[0][::-1]) if flipped else addp(plist[0])
    normal_one = obtain_normal3(*[vlist[i] for i in h_polygons[0][:3]])

    for idx in range(1, len(plist)):
        n = obtain_normal3(*[vlist[i] for i in plist[idx][:3]])
        if is_reasonably_opposite(n, normal_one):
            addp(plist[idx][::-1])
        else:
            addp(plist[idx])

    return h_polygons

def attracted(vlist, plist, flipped, attractors):
    '''
    - recalculate polygon direction based on attractors

    The median of the polygon is used to establish the location of a polygon.
    Each attractor is added to a kdtree, and for each median it's closest attractor
    is found and this relationship is used to determine if the polygon's direction
    should be flipped or not.

    '''

    #for each mesh:
    #    n = get normal of polygon, 
    #    if the average vector generated by connecting each vertex of the 
    #    polygon to the nearest attractor matches the found normal then keep index order, 
    #    else reverse
zeffii commented 7 years ago

some useful things scavenged from Blender's source, and restated in python.. GPL3 all of this

    def sub_v3_v3v3(a, b):
        return a[0]-b[0], a[1]-b[1], a[2]-b[2]

    def madd_v3_v3v3fl(a, b, f=1.0):
        return a[0]+b[0]*f, a[1]+b[1]*f, a[2]+b[2]*f

    def dot_v3v3(a, b):
        return a[0]*b[0]+a[1]*b[1]+a[2]*b[2]

    def isect_line_plane(l1, l2, plane_co, plane_no):

        u = l2[0]-l1[0], l2[1]-l1[1], l2[2]-l1[2]
        h = l1[0]-plane_co[0], l1[1]-plane_co[1], l1[2]-plane_co[2]
        dot = plane_no[0]*u[0] + plane_no[1]*u[1] + plane_no[2]*u[2]

        if abs(dot) > 1.0e-5:
            f = -(plane_no[0]*h[0] + plane_no[1]*h[1] + plane_no[2]*h[2]) / dot
            return l1[0]+u[0]*f, l1[1]+u[1]*f, l1[2]+u[2]*f        
        else:
            # The segment is parallel to plane
            return False
nortikin commented 7 years ago

mathutils?

zeffii commented 7 years ago

instead of casting to Vector. ( i don't know if this is going to be fast code...)

import math

def interp_v3_v3v3(a, b, t=0.5):
    if t == 0.0: return a
    elif t == 1.0: return b
    else:
        s = 1.0 - t
        return (s * a[0] + t * b[0], s * a[1] + t * b[1], s * a[2] + t * b[2])

def length(v):
    return math.sqrt((v[0] * v[0]) + (v[1] * v[1]) + (v[2] * v[2]))

def normalize(v):
    l = math.sqrt((v[0] * v[0]) + (v[1] * v[1]) + (v[2] * v[2]))
    return [v[0]/l, v[1]/l, v[2]/l]

def sub_v3_v3v3(a, b):
    return a[0]-b[0], a[1]-b[1], a[2]-b[2]

def madd_v3_v3v3fl(a, b, f=1.0):
    return a[0]+b[0]*f, a[1]+b[1]*f, a[2]+b[2]*f

def dot_v3v3(a, b):
    return a[0]*b[0]+a[1]*b[1]+a[2]*b[2]

def isect_line_plane(l1, l2, plane_co, plane_no):
    u = l2[0]-l1[0], l2[1]-l1[1], l2[2]-l1[2]
    h = l1[0]-plane_co[0], l1[1]-plane_co[1], l1[2]-plane_co[2]
    dot = plane_no[0]*u[0] + plane_no[1]*u[1] + plane_no[2]*u[2]

    if abs(dot) > 1.0e-5:
        f = -(plane_no[0]*h[0] + plane_no[1]*h[1] + plane_no[2]*h[2]) / dot
        return l1[0]+u[0]*f, l1[1]+u[1]*f, l1[2]+u[2]*f        
    else:
        # parallel to plane
        return False

def obtain_normal3(p1, p2, p3):
    return [
        ((p2[1]-p1[1])*(p3[2]-p1[2]))-((p2[2]-p1[2])*(p3[1]-p1[1])),
        ((p2[2]-p1[2])*(p3[0]-p1[0]))-((p2[0]-p1[0])*(p3[2]-p1[2])),
        ((p2[0]-p1[0])*(p3[1]-p1[1]))-((p2[1]-p1[1])*(p3[0]-p1[0]))
    ]

def mean(verts):
    vx, vy, vz = 0, 0, 0
    for v in verts:
        vx += v[0]
        vy += v[1]
        vz += v[2]
    numverts = len(verts) 
    return vx/numverts, vy/numverts, vz/numverts

def perform_recalcs(input_data):
    num_data_channels = len(input_data)
    func = {3: homogenous, 4: attracted}.get(num_data_channels)
    return func(*input_data) if func else []

def is_reasonably_opposite(n, normal_one):
    return dot_v3v3(normalized(n), normalized(normal_one)) < 0.0

def homogenous(vlist, plist, flipped):
    '''
    - recalculate homogenous

    assume all normals are perpendicular to a single plane, but some normals 
    are flipped. This will turn them all to face the direction of the first found
    normal (or reverse if flipped is True).

    useful for islands of potentially flipped normals, where direction can not be 
    usefully inferred by vacinity.

    unfortunately, this doesn't work great on concave polygons..
    '''

    h_polygons = []
    addp = h_polygons.append

    addp(plist[0][::-1]) if flipped else addp(plist[0])
    normal_one = obtain_normal3(*[vlist[i] for i in h_polygons[0][:3]])

    for idx in range(1, len(plist)):
        n = obtain_normal3(*[vlist[i] for i in plist[idx][:3]])
        if is_reasonably_opposite(n, normal_one):
            addp(plist[idx][::-1])
        else:
            addp(plist[idx])

    return h_polygons

def attracted(vlist, plist, flipped, attractors):
    '''
    - recalculate polygon direction based on attractors

    The median of the polygon is used to establish the location of a polygon.
    Each attractor is added to a kdtree, and for each median it's closest attractor
    is found and this relationship is used to determine if the polygon's direction
    should be flipped or not.

    use -----------

    for polygon of plist:
        n = get normal of polygon, 
        if the average vector generated by connecting each vertex of the 
        polygon to the nearest attractor matches the found normal then keep index order, 
        else reverse
    '''

    h_polygons = []
    failures = []

    if not all(plist, attractors):
        return [], []

    size = len(vlist)
    kd = mathutils.kdtree.KDTree(size)

    for i, point in enumerate(attractors):
        kd.insert(point, i)
    kd.balance()

    for polygon in plist:
        avg = mean([vlist[vidx] for vidx in polygon])
        plane_no = obtain_normal3(*[vlist[vidx] for vidx in polygon[:3]])
        co, index, dist = kd.find_n(avg, 1)[0]
        p = isect_line_plane(l1=avg, l2=co, avg=avg, plane_no=plane_no)
        if not p:
            failures.append(polygon)
        else:
            n  = normalized(sub_v3_v3v3(plane_no, avg))
            b = is_reasonably_opposite(n, co)
            h_polygons.append(polygon if b else polygon[::-1])

    return h_polygons, failures
zeffii commented 7 years ago

can make things inline :) just want to test it first..

zeffii commented 7 years ago

code was finished some time last week, didn't get a chance to make a node out of it. close for now..