Pyro 2 Work Flow Examples

I totally forgot I had these laying around.. Thought I’d share them.
Simple examples on how to use some of the Pyro 2 functionality, such as inheriting object motion, object collision, clustering etc.
Not very visual but could be useful. Small post-it notes explain the various subsets present in the netwerk.

Note that the examples only work with Houdini 12 +.

Get the package HERE

Automated Edge Flip Operator

Last week I ran in to a long lasting Houdini frustration, the lack of proper Edge support.
We get the control over various primitive types, points and verts, but not Edges.

This becomes frustrating when dealing with complex primitive sorting and optimization techniques. For that purpose I decided to create my own Edge object, a rather simple wrapper that let’s me create Edges from points, verts and primitives inside Houdini. An edge can be constructed using two points, a primitive or 2 verts. Functionality includes comparing angles, getting the length, comparing other existing edges on point presence, not order etc.
Just for convenience 🙂

One of the first things I needed it for was an automatic edge flip operator. The SOP samples primitive data and flips edges based on the optimal angle of neighbouring prims. Similar to the Delauney edge orientation principle. The flipping of an edge is only possible when the union of the two incident faces forms a convex quadrilateral. Flip Constraints are based on the curvature of a reference mesh, optimizing the triangle edge orientation on harsh edges or cutoff points.

To determine if the edge needs to be flipped, both edge lengths are computed, together with the maximum shared angle. If the angle (thus the edge length) is higher for the incident face, the edge can be flipped. Otherwise the angle is compared based on the maximum difference between the point normals and prim to point normals. If the difference comes close to 0, the edge follows curvature. If it is higher then the specified threshold, the edge should be flipped. (conforming the proposed mesh orientation and curvature).

One of the implementation issues: how do you keep track of the mesh when you flip an edge? Flipping an edge makes the next edge (sharing one of the previous edge primitives) most likely invalid, or a possible candidate. You could solve this by recursively update the mesh as you flip an edge, keeping track of the deforming state and marking areas that have been computed. This would be rather slow and tedious to incorporate. I ended up marking the edge that should be flipped, any new edge being sampled is being checked against any shared edge points. If one of the shared edge points match a flipped edge, it isn’t considered. Running the operator a couple of times would asure most edges will align. Most likely ater two computational rounds. I also added a “dirty” mode, not taking into account previously flipped edges. (works surprisingly well)

The actual operator doesn’t flip the edges. Houdini’s Edge Flip Sop does that for me. I simply count the total number of edges to be flipped and save the Edge Points as a Detail String Attribute.
The count is used to specify the number of iterations the loop should make and inside I simply sample each saved edge. It is a limitation I had to work with in Houdini.

A handy tip on speed and flexibility might be to pre-sample your data-set in to clusters. Clusters are formed based on any Vector3 attribute and are weighted over X recursions. Finding the most optimal distribution as it is being computed. Creating clusters (based on Position for example) speeds up the computing of the edges to be flipped, because the scope is limited to only a small island of primitives. In comparison to the entire mesh. Similar to using a Kd-Tree, but simpler. Read up on the algorithm here and here

Now for some code, the edge class should be easy to recreate based on what was told and can be seen in the following snippet:

from Triangulate.hedge import Edge


# Helps me debug my code
debug = hou.pwd().evalParm("debug_mode")
export_name = hou.pwd().evalParm("edge_name")
flip_count_name = hou.pwd().evalParm("count_name")
export_groups = hou.pwd().evalParm("export_groups")
compare_angle = hou.pwd().evalParm("compare_angle")
angle_attrib_name = hou.pwd().evalParm("angle_attribute")
constrain_attrib_name = hou.pwd().evalParm("constrain_attribute")
constrain_angle = hou.pwd().evalParm("constrain_angle")
min_angle = hou.pwd().evalParm("min_angle")
max_angle = hou.pwd().evalParm("max_angle")
constrain_cut = hou.pwd().evalParm("constrain_cutoff")
safe_mode = hou.pwd().evalParm("safe_mode")

# This code is called when instances of this SOP cook.
node = hou.pwd()
geo = node.geometry()
prims = geo.prims()
points = geo.points()
edge_list = []
shared_edges = []
flipped_edges = []
flipped_points = []

Necessary Attributes---------------------------------------------------------------------------------------------------

# Find the count attribute or create it otherwise
count_attrib = geo.findGlobalAttrib(flip_count_name)
if count_attrib is None:
    count_attrib = geo.addAttrib(hou.attribType.Global, flip_count_name, 0)
    geo.setGlobalAttribValue(count_attrib, 0)

angle_attrib = geo.findPointAttrib(angle_attrib_name)
if angle_attrib is None and compare_angle:
    raise hou.NodeError("Can't find angle comparison point attribute: %s" % angle_attrib_name)

constrain_attrib = geo.findPointAttrib(constrain_attrib_name)
if constrain_attrib is None and constrain_angle and compare_angle:
    raise hou.NodeError("Can't find constrain attribute: %s" % constrain_attrib_name)


def Console(inString):
    Helps me to debug my code

    if debug:
        print inString

def GetPrimPoints(inPrim):
    Get's the points out of a primitive

    point_list = []
    for vert in inPrim.vertices():
    return point_list

def GetHighestAngle(inNormals):
    Compares vectors to get the highest angle

    max_angle = 0
    for n in xrange(len(inNormals)):
        n1 = inNormals[n]
        for ni in range(n+1, len(inNormals)):
            n2 = inNormals[ni]
            angle = n1.angleTo(n2)
            if angle > max_angle:
                max_angle = angle

    return max_angle

def GetHighestPrimAngle(inNormals, inPrimNormal):
    Compares the normals against the primitive

    max_angle = 0
    for normal in inNormals:
        angle = normal.angleTo(inPrimNormal)
        if angle > max_angle:
            max_angle = angle

    return max_angle

def FindUnsharedPoints(inPrimOne, inPrimTwo, inEdge):
    Returns a list of points that isn't shared between the primitives

    unshared_points = []
    prim_one_points = GetPrimPoints(inPrimOne)
    prim_two_points = GetPrimPoints(inPrimTwo)

    # compare first primitve
    for point in prim_one_points:
        if point in inEdge.points:

    # compare second primitive
    for point in prim_two_points:
        if point in inEdge.points:

    return unshared_points

def SampleEdges(inPrim):
    Computes the primitive edges
    New edges are appended to the edge_list

    Console("sampling prim: %d" % inPrim.number())

    #Fetch the vertices
    point_list = GetPrimPoints(inPrim)

    # Now calculate the edge
    amp = len(point_list)
    for i in range(amp):
        if i < (amp-1):
            point_one = point_list[i]
            point_two = point_list[i+1]
            point_one = point_list[-1:][0]
            point_two = point_list[0]

        edge = Edge(point_one, point_two)

        # if the edge is part of already existing edges, it shared a primitive
        # i know this is dirty but i'd rather have this then do 2 lookups :(
            Console("shared edge: %s" % edge)
            existing_edge = edge_list[edge_list.index(edge)]
        except ValueError:
            Console("new edge: %s" % edge)

def CompareAngle(inPrimOne, inPrimTwo, inEdge, inUnsharedPoints):
    Compares the angle of the point normal angle against primitive normal
    Returns true if the edge needs to be inverted...
    The comparison is based on the maximum difference between the prim point normals,
    and the prim point normals compared with the prim normal. If the difference is rather
    large on both sides, consider the edge to be flipped

    An extra constrain attribute can be specified. Works best with curvature.
    If the curvature of the unshared points is higher then the combined edge points,
    It is allowed to flip, together with the angle constrain

    prim_one_points = GetPrimPoints(inPrimOne)
    prim_two_points = GetPrimPoints(inPrimTwo)

    # get the point normals
    prim_one_normals = []
    prim_two_normals = []
    for point in prim_one_points:
    for point in prim_two_points:

    # figure out the angles
    max_angle_one = GetHighestAngle(prim_one_normals)                                  # the highest measured angle between the point normals
    max_angle_two = GetHighestAngle(prim_two_normals)                                
    max_prim_angle_one = GetHighestPrimAngle(prim_one_normals, inPrimOne.normal())     # the highest measued angle between points and prim normal
    max_prim_angle_two = GetHighestPrimAngle(prim_two_normals, inPrimTwo.normal())

    angle_diff_one = abs(max_angle_one - max_prim_angle_one)                           # difference between prim 1 normals      
    angle_diff_two = abs(max_angle_two - max_prim_angle_two)                           # difference between prim 2 normals

    max_value = max(angle_diff_one, angle_diff_two)                                    # get the higest of the two                            
    min_value = min(angle_diff_one, angle_diff_two)                                    # get the lowest of the two
    Console("normal angle one: %s" % max_angle_one)                                    # some debug stuff
    Console("normal angle two: %s" % max_angle_two)
    Console("prim   angle one: %s" % max_prim_angle_one)
    Console("prim   angle two: %s" % max_prim_angle_two)
    Console("angle diff one: %s" % angle_diff_one)
    Console("angle diff two: %s" % angle_diff_two)

    # figure out the additional attribute
    if constrain_angle:
        com_original = inEdge.pointOne.attribValue(constrain_attrib) + inEdge.pointTwo.attribValue(constrain_attrib)            # original edge point attribute values
        com_diff = inUnsharedPoints[0].attribValue(constrain_attrib) + inUnsharedPoints[1].attribValue(constrain_attrib)        # possible new edge attribute values
        constrain = com_diff / 2.0                                                                                              # average to be compared against attribute cutoff

        if com_diff > com_original and constrain > constrain_cut and max_value >= max_angle and min_value >= min_angle:
            return True
        if max_value >= max_angle and min_value >= min_angle:
            return True

    return False

def FlipEdge(inEdge):
    Checks if the edge needs to be flipped or not

    if (safe_mode) and (inEdge.pointOne in flipped_points or inEdge.pointTwo in flipped_points):
        Console("ALREADY IN FLIPPED EDGES: %s" % inEdge)
        return -1

    prim_one = prims[inEdge.edgePrims[0]]
    prim_two = prims[inEdge.edgePrims[1]]
    # for both primitives find the point that isn't shared and figure out if it should be flipped
    unshared_points = FindUnsharedPoints(prim_one, prim_two, inEdge)
    if len(unshared_points) != 2:
        ("%s has more more or less then 2 unshared points" % inEdge)
        return -1

    Console("unshared points: %s, %s" % (unshared_points[0].number(), unshared_points[1].number()))
    # check edge length
    edge_two = Edge(unshared_points[0], unshared_points[1])
    edgelength_one = inEdge.edgeLength
    edgelength_two = edge_two.edgeLength
    if edgelength_two < edgelength_one:
        Console("FLIPPING EDGE!")
        if safe_mode:
        return 1

    # check angle
    elif (compare_angle) and (CompareAngle(prim_one, prim_two, inEdge, unshared_points)):
        Console("FLIPPING EDGE!")
        if safe_mode:
        return 1

def Finalize():
    Creates the detail attributes and groups (if chosen)
    The counter specifies the amount of edges to be flipped

    count = 0
    for edge in flipped_edges:

        # set or create the detail attribute
        edge_attrib_name = "%s_%d" % (export_name, count)
        edge_attrib = geo.findGlobalAttrib(edge_attrib_name)
        if edge_attrib is None:
            edge_attrib = geo.addAttrib(hou.attribType.Global, edge_attrib_name, "none")
        geo.setGlobalAttribValue(edge_attrib,"p%d-%d" % (edge.pointNumberOne, edge.pointNumberTwo))
        # see if we export the group
        if export_groups:
            group = geo.createPointGroup("%s%d_%d" % (export_name, edge.pointNumberOne, edge.pointNumberTwo))
            for point in edge.points:

        count += 1

    # set the global attribute value
    geo.setGlobalAttribValue(count_attrib, count)


# For every prim construct 3 edges
# See if the edges have multiple primitives shared
# If multiple primitives are shared, do the comparison

# Calculate shared edges
for prim in prims:

# See if they need to be flipped
for edge in shared_edges:
    Console("sampling edge: %s" % edge)

# Set up the edges to be flipped

Primitive Normal Orientation

Ever wondered what “Orient Polygons” does? It tells me that it winds all the polygons in the same ‘direction’. What direction is that?
And how does it know (determine) what polygons to reverse or flip. I find that desription utterly useless and needed more control over my prim normals in Houdini (defining direction using a set of shared points and their surface orientaiton) for a project i’m working on.

Based on those requirementrs I made a simple orient prim normal operator that uses a vector point attribute to test the primitive orientation. If the op determines the prim normal needs to be flipped, it’s added to a primitive normal group that can be reversed using the reverse SOP.

import math

# This code is called when instances of this SOP cook.
node = hou.pwd()
geo = node.geometry()

# Evaluate the parameters
prim_group = hou.pwd().evalParm("primitive_group")          #primitive group to sample from
normal_name = hou.pwd().evalParm("sample_attribute")        #vector attribute used for averaging and comparison
export_group = hou.pwd().evalParm("export_group")           #primitive group wrong faced normals get exported to
angle_threshold = hou.pwd().evalParm("angle_offset")     #maximum allowed angle

# Finds the primitive group
def findPrimGroup(inName):

    for group in geo.primGroups():
        if == inName:
            return group.prims()

    return None

Sample Information------------------------------------------------------------------------------------------

# Fetch the primitives and sample input information
if len(prim_group) is not 0:
    prims = findPrimGroup(prim_group)
    if prims is None:
        raise hou.NodeError("Can't find primitive group: %s" % prim_group)

    prims = geo.prims()

if len(prims) < 1:
    raise hou.NodeError("Not enough primitives provided! Need at least one..")

if len(export_group) < 1:
    raise hou.NodeError("No export group found")

if export_group in geo.primGroups():
    raise hou.NodeError("Primitive group '%s' already exists" % export_group)

export_group = geo.createPrimGroup(export_group)

normal_attrib = geo.findPointAttrib(normal_name)
if normal_attrib is None:
    raise hou.NodeError("Lookup attribute '%s' does not exist" % normal_name)


# Now sample the normals
for prim in prims:
    # get the prim verts
    verts = prim.vertices()

    # from every vert get the associated point normals
    # note that every point normal is normalized
    prim_point_normals = []
    for v in verts:
        vector = hou.Vector3(v.point().attribValue(normal_attrib))

    # accumulate point normals
    accum = hou.Vector3(0,0,0)
    for n in prim_point_normals:
        for i in range(3):

    # get the average
    count = len(prim_point_normals)
    for i in xrange(len(accum)):
        accum[i] /= count

    # get primitive normal
    prim_normal = hou.Vector3(prim.normal())
    inv_prim_normal = -1*prim_normal

    # get the angle
    current = accum.angleTo(prim_normal)
    inversed = accum.angleTo(inv_prim_normal)

    if inversed < (current+angle_threshold):

Houdini Pyro2 Dev Tests

While working on Pyro2 at Side Effects, some time was spend on simming and rendering dev tests. After incorporating a new feature (usability or tech related), a test drive was necessary.

Houdini 12’s Demo Video already shows what is capable using the new set of tools. But I though it might be interesting to see what was rendered trying to improve the tools.

All small tests were rendered using Mantra and simmed using the tools being worked on at that time. The Pyro 1 shader was used, as version 2 wasn’t available.

Note that leaving Velocity Blur off, some flames might appear a bit ‘blobsy’. The streaking effect will work better turning velocity blur on.

Render settings:

  • Micro Polygon Rendering
  • Volume Filter Width, 1 or 0.75
  • 1 or 2k Transparent Depth Map shadows
  • Max of 3 spot lights, classic 3 way lightning set-up
  • Max 15 minutes a frame, 720×480
  • In case of Clustering, Delayed Load archives were used (more memory efficient)
  • 3-6 pixel samples
  • Gaussian 2×2 Pixel Filter
  • Temperature was used as a transparency mask for Heat
  • No Noise
  • No Velocity Blur

Burning Sphere

Burning Grid

Heavy Fire

Moving Fire Ball

Clustered Tables

Small Smoke Trail

Explosion Playblast (viewport)

Point Cloud Meshing HOM

I’ve been working on a point cloud meshing method, incorporating this Paper

The algorithm has been extended using a custom nearest neighbour and Forward Facing Normal Flow operator written in VEX. I might be able to post those soon.

Anyway, the surface reconstruction is handled by a custom Python Operator, connecting neighbours, sampling edges, defining edge normals etc.

The most optimal distribution is searched for, empty spots are recursively filled. This is a work in progress and I won’t be able to post the complete functioning operator. But this snippet of code might help (tip of the iceberg ;)).

import Triangulate.hedge
import Triangulate.hpoly


from Triangulate.hedge import Edge
from Triangulate.hpoly import TriPoly
import sys
import math

# This code is called when instances of this SOP cook.
node = hou.pwd()
geo = node.geometry()


# Set up some globals
neigh_name = hou.pwd().evalParm("neighbour_name")                   # neighbour attribute name
angle_name = hou.pwd().evalParm("angle_name")                       # angle attribute name
distance_name = hou.pwd().evalParm("distance_name")                 # distance attribute name

edge_group_name = hou.pwd().evalParm("edge_group_name")             # group of initial edge points
neigh_amount = hou.pwd().evalParm("neighbours")                     # amount of neighbours
max_iter = hou.pwd().evalParm("max_iterations")                     # maximum amount of iterations
point_prior = hou.pwd().parm("point_priorization").evalAsString()   # point sampling method
max_connections = hou.pwd().evalParm("max_connections")             # max allowed point connections
conn_constrain = hou.pwd().parm("constrain_method").eval()          # point constrain method
double_constrain = hou.pwd().parm("double_reduction").eval()        # double prims constrain method
orient_polygons = hou.pwd().parm("orient_polygons").eval()          # if polygons need to be oriented
predefined_edges = hou.pwd().evalParm("predefined")                 # if edges are predefined in groups

num_connections_name = hou.pwd().evalParm("connections_attrib")     # name of number of connections
num_associations_name = hou.pwd().evalParm("association_attrib")    # name of the association attribute
unresolved_name = hou.pwd().evalParm("unresolved_attribute")        # name of the unresolved edge attribute
orient_name = hou.pwd().evalParm("orient_attribute")                # name of the orientation attribute

orient_group_name = hou.pwd().evalParm("orient_export_name")        # name of the orient polygons group
unresolved_group_name = hou.pwd().evalParm("unresolved_export_name")# name of the unresolved prims group
unresolved_edges_name = hou.pwd().evalParm("unresolved_edges_name")    

DEBUG = hou.pwd().evalParm("debug")

neigh_attributes = []                                               # list of neighbour point attributes
connected_points = []                                               # list of point numbers that have been used!
active_edges = []                                                   # holds the currently active edges
unresolved_edges = []                                               # holds the currently unresolved edges
boundary_edges = []                                                 # holds if an edge is on a boundary
bound_polygons = []                                                 # list of polygons created

geo = hou.pwd().geometry()                                          # reference to incoming geo
points = geo.points()                                               # original points


# Add the number of connections attribute
num_connections_attrib = geo.findPointAttrib(num_connections_name)
if num_connections_attrib is None:
    num_connections_attrib = geo.addAttrib(hou.attribType.Point, num_connections_name, 0)

# Add the number of associations attribute
num_associations_attrib = geo.findPointAttrib(num_associations_name)
if num_associations_attrib is None:
    num_associations_attrib = geo.addAttrib(hou.attribType.Point, num_associations_name, 0)

# Add the point resolved attribute
resolved_attrib = geo.findPointAttrib(unresolved_name)
if resolved_attrib is None:
    resolved_attrib = geo.addAttrib(hou.attribType.Point, unresolved_name, 0)

resolved_prim_attrib = geo.findPrimAttrib(unresolved_name)
if resolved_prim_attrib is None:
    resolved_prim_attrib = geo.addAttrib(hou.attribType.Prim, unresolved_name,0)


#special debug method
def Console(inString):
    Helper function showing debug information if switched on

    if DEBUG:
        print inString
def GetPoint(inNumber):
    Returns a point by index

    return points[inNumber]

def GetNumberOfAssociations(inPoint):
    Returns the number of point associations

    return inPoint.numberOfAssociations

def GetNumberOfConnections(inPoint):
    Returns the number of point connections

    return inPoint.numberOfConnections

def PolyAssociatedPresence(inEdge, inPoint):
    Checks by heritage if the poly is already created

    if inEdge.heritagePoint == inPoint:
        return True
    return False

def PolyPresence(inEdge, inPoint):
    Check in the bound polygon list if the poly was already created

    if TriPoly([inEdge.pointOne, inEdge.pointTwo, inPoint]) in bound_polygons:
        return True
    return False

def findPrimGroup(inName):
    Finds a primitive group with a specific name
    Returns none if not found!

    for group in geo.primGroups():
        if == inName:
            return group.prims()

    return None

# get's the edge group points
def GetSeedPoints(inGroupName):
    Returns the points of the specified point group

    return_group = None
    for group in geo.pointGroups():
        if == inGroupName:
            if len(group.points()) > 2:
                pgroup = []
                for point in group.points():                            # Search for the point in the original points and return those!
                return pgroup
                raise hou.NodeError( "Not enough points in edge group, minimum of two required!")
                return -1

    raise hou.NodeError( "Initial point group not found or not enough points: %s" % inGroupName)
    return -1

# adds some sampling attributes
def AddPointAttributes(inPoints):
    Adds number of connections and isEdge attribute to points

    for point in inPoints:
        point.numberOfConnections = 0                      #set number of connections to 0
        point.numberOfAssociations = 0                     #set number of point associations to 0
        point.isEdgePoint = False                          #see if it's an edge point or not
        #point.isUnresolved = 1                             #set's the point to be unresolved (share a min of 2 edges)

# creates seed edges
def CreateSeedEdges(inPoints):
    Based on the seed points, create some edges to start with
    These edges are directly set to be active and to be an edge

    for pnum in range(len(inPoints)):
        if pnum <= len(inPoints)-2:
            edge = Edge(inPoints[pnum], inPoints[pnum+1])
            edge = Edge(inPoints[pnum], inPoints[0])

        active_edges.append(edge)                           #add to active edge
        boundary_edges.append(edge)                         #add to bound  edge
        connected_points.append(inPoints[pnum].number())    #add as connected point
        inPoints[pnum].isEdgePoint = True                   #set it to be an edge point
        inPoints[pnum].numberOfConnections = 2              #set it to already have 2 connections
        inPoints[pnum].numberOfAssociations = 0             #set number of point associations to 0
def CreatePredefinedEdges(inGroupName):
    Creates a list of predefined edges that were submitted in groups
    Based on the initial group name, appendices are sampled and added as an edge

    for group in geo.pointGroups():
        if inGroupName in
            # check to see if we're dealing with 2 points
            group_points = group.points()
            if len(group_points) != 2 :
                raise hou.NodeError("Error in group: %s Only predefined groups of two points can be sampled" % group)
            edge = Edge(group_points[0], group_points[1])

# initializes the neigbours
def InitializeNeighbourAttributes(inGeo, inName):
    initialized a list of neighbour attributes

    rlist = []
    for n in range(neigh_amount):
        attr = inGeo.findPointAttrib("%s%d" % (inName,(n+1)))
        if attr != None:
            print "can't find neighbour attribute: %s%d" % (inName, (n+1))
    return rlist

# initialize export group and attribute
def InitializeGroup(inName):
    Initializes the export group and orientation attribute

    if len(inName) < 1:
        raise hou.NodeError("Please specify an export group, no name found")

    if inName in geo.primGroups():
        raise hou.NodeError("Primitive group '%s' already exists" % inName)

    return geo.createPrimGroup(inName)

# initialize orientation attribute
def InitializeOrientAttrib(inName):
    Initializes the poly orientation attribute

    orient_attrib = geo.findPointAttrib(inName)
    if orient_polygons and orient_attrib is None:
        raise hou.NodeError("Can't find orientation attribute: %s" % orient_name)
    return orient_attrib

# finalizes the attributes
def Finalize(inUnresolvedGroup):
    Transfers point properties on to attributes

    # Set the number of connections attribute
    for point in points:
        point.setAttribValue(num_connections_attrib, point.numberOfConnections)
        point.setAttribValue(num_associations_attrib, point.numberOfAssociations)

    edge_group_counter = 0

    # Append unresolved primitives to unresolved group and set attribute
    for edge in unresolved_edges:
        Console( "Unresolved Edge: %s, Prim: %s" % (edge, edge.edgePrims))
        # Create the edge group
        group = geo.createPointGroup("%s_%d" % (unresolved_edges_name, edge_group_counter))
        edge_group_counter += 1
        for prim in edge.edgePrims:
            prim = geo.iterPrims()[prim]
            prim.setAttribValue(resolved_prim_attrib, 1)
        for point in edge.points:
            point.setAttribValue(resolved_attrib, 1)
    # Append unresolved active primitives to unresolved group and set attribute
    for edge in active_edges:
        Console("Unresolved Active Edge: %s, Prim: %s" % (edge, edge.edgePrims))
        group = geo.createPointGroup("%s_%d" % (unresolved_edges_name, edge_group_counter))
        edge_group_counter += 1
        for prim in edge.edgePrims:
            prim = geo.iterPrims()[prim]
            prim.setAttribValue(resolved_prim_attrib, 1)
        for point in edge.points:
            point.setAttribValue(resolved_attrib, 1)

# updates the edges
def UpdateActiveEdges(inEdges):
    Verifies if the incoming edges are already being used
    If used, remove from the active edge list, as the poly is most likely closed
    Also increment number of connections

    for edge in inEdges:
        if edge in active_edges:                                       #if in active_edges, remove
            Console("edge was in active edges: %s" % str(edge))

        if edge in unresolved_edges:                                   #if in unresolved edges, remove
            Console("edge was in unresolved edges: %s" % str(edge))
        active_edges.append(edge)                                      #otherwise increment connections and allow for sampling
        Console("Edge NOT in active edges: %s" % str(edge))
        for point in edge.points:
            point.numberOfConnections += 1

# update point associations
def UpdateAssociations(inPoints):
    increments the 3 vertices of the just created triangle

    for point in inPoints:
        point.numberOfAssociations += 1

# orient polygons
def OrientPolygons(inPoints, inPoly):
    compares the point normals agains prim normal

    prim_point_normals = []
    for point in inPoints:
        vector = hou.Vector3(point.attribValue(orient_attrib))
    # combine normals    
    accum = hou.Vector3(0,0,0)
    for n in prim_point_normals:
        for i in range(3):
    # get the average
    count = len(prim_point_normals)
    for i in xrange(len(accum)):
        accum[i] /= count
    # get the prim normal
    prim_normal = hou.Vector3(inPoly.normal())
    inv_prim_normal = -1*prim_normal
    # compare
    if accum.angleTo(inv_prim_normal) < accum.angleTo(prim_normal):

# sort polygon verts        
def SortVertices(inPoints):
    Sorts the incoming points clockwise

    order_array = [None,None,None]          # start off with an empty array          
    xmin = sys.float_info.max
    zmin = sys.float_info.max
    for point in inPoints:                  # find first point
        posx = point.position()[0]
        if posx < xmin:
            xmin = posx
            order_array[0] = point
    for point in inPoints:                  # find second point
        posz = point.position()[2]
        if posz < zmin:
            zmin = posz
            order_array[1] = point
    inPoints.remove(order_array[1])         # append third point
    order_array[2] = inPoints[0]
    return order_array

# Create a polygon out of a previous edge and new point
def CreatePolygon(inEdge, inPoint):
    Create A poly out of an edge and a point
    TODO:  Add a edge normal to every edge created, + point if necess ary

    poly = geo.createPolygon()
    if inEdge.pointOne.isEdgePoint and inEdge.pointTwo.isEdgePoint:
        order_array = [inEdge.pointOne, inEdge.pointTwo, inPoint]
        order_array = SortVertices([inEdge.pointOne, inEdge.pointTwo, inPoint])
    # add the poly as a TriPoly to our list of polygons
    # check poly orientations
    if orient_polygons:
        OrientPolygons(order_array, poly)

    # create the new edges
    edge_list = [Edge(inEdge.pointOne, inPoint, inEdge.pointTwo), Edge(inEdge.pointTwo, inPoint, inEdge.pointOne), inEdge]
    # The edge now needs it's associated polygon
    for edge in edge_list:

    # see if the edge already exists, if so, no active edges should be added
    # increase number of references to point
    UpdateAssociations([inEdge.pointOne, inEdge.pointTwo, inPoint])

def FilterSharedNeighbour(inDictOne, inDictTwo, inEdge):
    iterates through the points in the dictionary, searching for the best shared neighbour.
    sampling has two options: `minimum` or `closest`. Minimum = abs(x-y), closest = x+y
    returns None if no suitable point is found
    TODO: Add pointcheck that point is in between other points

    closest_point = None
    current_dist = sys.float_info.max
    conn_dict = {0: GetNumberOfConnections, 1: GetNumberOfAssociations}                         # holds the two connection comparison methods
    double_dict = {0: PolyAssociatedPresence, 1: PolyPresence}

    for key_point in inDictOne:
        if inDictTwo.has_key(key_point):
            sample_point = GetPoint(key_point)                                                  # returns the associated point

            connections = conn_dict[conn_constrain](sample_point)

            if sample_point.isEdgePoint or connections >= max_connections:                      # compares the neighbour based on edge, connectivity
                Console("invalid point: %s" % key_point)

            if double_dict[double_constrain](inEdge, sample_point):                             # compares the neighbour based on poly existance
                Console("creates similar poly: %s" % key_point)

            p2_dist = inDictTwo[key_point][0]                                                   # get the point two distance
            p1_dist = inDictOne[key_point][0]                                                   # get the point one distance
            if point_prior == "minimal":                                                        # choose sampling method
                dist = abs(p1_dist - p2_dist)
                dist = p1_dist + p2_dist
            Console("shared point: %s, distance: %s" % (key_point, dist))

            if dist < current_dist:
                closest_point = sample_point
                current_dist = dist
    return closest_point

def GetEdgeNeighbour(inEdge):
    Finds the most suitable edge point to create a triangle with
    In case multiple points are found, the one with the lowest distributed distance
    is chosen.

    # sample the values
    nlist_one = {}
    nlist_two = {}

    point_one = inEdge.pointOne
    point_two = inEdge.pointTwo

    Console("sampling edge: %s" % str(inEdge))

    # we want to construct a dict with: `pointnum:[dist, angle]`
    counter = 0
    for attribute in neigh_attributes:
        # first get the point numbers as keys
        pnum_one = point_one.attribValue(attribute)
        pnum_two = point_two.attribValue(attribute)
        # sample the distance and angle for every neighbour point
        pdist_one = point_one.attribValue(dist_attributes[counter])
        pdist_two = point_two.attribValue(dist_attributes[counter])
        pangle_one = point_one.attribValue(angle_attributes[counter])
        pangle_two = point_two.attribValue(angle_attributes[counter])
        nlist_one[pnum_one] = [pdist_one, pangle_one]
        nlist_two[pnum_two] = [pdist_two, pangle_two]

        counter += 1
    # now we have the point numbers + distance / angle, do the comparison
    closest_point = FilterSharedNeighbour(nlist_one, nlist_two, inEdge)

    if closest_point is not None:
        Console("closest shared point: %s" % closest_point.number())

    return closest_point


# Initialize the unresolved poly group
unresolved_group = InitializeGroup(unresolved_group_name)

# If an orientation is specified, we want sample the export group
if orient_polygons:
    export_group = InitializeGroup(orient_group_name)
    orient_attrib = InitializeOrientAttrib(orient_name)

# Initialize the points with the necessary attributes

# Initialize the neighbour attributes used for sampling best point to angle fit
neigh_attributes = InitializeNeighbourAttributes(geo, neigh_name)
dist_attributes = InitializeNeighbourAttributes(geo, distance_name)
angle_attributes = InitializeNeighbourAttributes(geo, angle_name)

# Now create a start list of active edges
if not predefined_edges:

# Iterate over the active edges and connect
count = 0
while len(active_edges) > 0 and count < max_iter:

    current_edge = active_edges[0]
    closest_point = GetEdgeNeighbour(current_edge)
    if closest_point != None:
        CreatePolygon(current_edge, closest_point)
    count += 1



Houdini Edge Normals

Currently I’m working on various point cloud meshing methods.

For these methods to work I needed to calculate perpendicular edge normals.

Quick python op for calculating specific edge normals.. When applied creates a triangle + edge normals

# This code is called when instances of this SOP cook.
node = hou.pwd()
geo = node.geometry()

# calculates the edge normal using the prim normal
def calculateEdge(inPoint1, inPoint2, inPolyNormal):

    diff = inPoint2.position() - inPoint1.position()

    edge_point_position = (diff*0.5) + inPoint1.position()
    edge_point = geo.createPoint()

    diff_normalized = diff.normalized()
    poly_normal = inPolyNormal.normalized()
    point_normal = poly_normal.cross(diff_normalized)

    edge_point.setAttribValue(attrib, point_normal)
    edge_point.setAttribValue(is_edge_attrib, 1)

# Evaluate offset parameter
pp1 = hou.pwd().evalParmTuple("point_one")
pp2 = hou.pwd().evalParmTuple("point_two")
pp3 = hou.pwd().evalParmTuple("point_three")

# Add the normal attribute
attrib = geo.addAttrib(hou.attribType.Point, "edge_normal", (0.0,0.0,0.0), transform_as_normal=True)
is_edge_attrib = geo.addAttrib(hou.attribType.Point, "edge_point", 0)

# Create the points
p1 = geo.createPoint()
p2 = geo.createPoint()
p3 = geo.createPoint()


point_list = [p1,p2,p3]

# Create Tri poly
poly = geo.createPolygon()

# Find the normal for ever edge
for i in range(len(point_list)):

    if i < (len(point_list)-1):
        point_one = point_list[i]
        point_two = point_list[i+1]
        point_one = point_list[-1:][0]
        point_two = point_list[0]

    calculateEdge(point_one, point_two, poly.normal())

Gaussian Blur Filter C++

I got asked to make some new blur filters.

The available convolution filters turned out to be rather slow and a set of new ones was requested. Getting to know the specific PDK (plugin development kit) was tricky, writing the plug-ins on the other end was a lot of fun. Below there is a snippet of code on how to write a Gaussian and Box blur kernel in C++.

Note that the height field input parameter (HField *inHeightMap and BuildContext &inContext)   can be replaced with any other (pixel) matrix using for example Devil or FreeImage.

The BuildContext is used to generate a second image and is used as a temporary pixel placeholder.

Now for the code:

I’ll start by defining the Gaussian Kernel. When called this kernel function returns an array of floats that define the actual per pixel scalar values.

// Calculates a 1d gaussian bell shaped kernel
float* GBlur::ComputeGaussianKernel(const int inRadius, const float inWeight)
    int mem_amount = (inRadius*2)+1;
    float* gaussian_kernel = (float*)malloc(mem_amount*sizeof(float));

    float twoRadiusSquaredRecip = 1.0 / (2.0 * inRadius * inRadius);
    float sqrtTwoPiTimesRadiusRecip = 1.0 / (sqrt(2.0 * PI) * inRadius);
    float radiusModifier = inWeight;

    // Create Gaussian Kernel
    int r = -inRadius;
    float sum = 0.0f;
    for (int i = 0; i < mem_amount; i++)
        float x = r * radiusModifier;
        x *= x;
        float v = sqrtTwoPiTimesRadiusRecip * exp(-x * twoRadiusSquaredRecip);
        gaussian_kernel[i] = v;

    // Normalize distribution
    float div = sum;
    for (int i = 0; i < mem_amount; i++)
        gaussian_kernel[i] /= div;

    return gaussian_kernel;

Now this kernel can be used to blur an image. The image comes in as a an array of floats.
This array represents a greyscale image. A colored image would have the RGB(A) components blurred.

// Calculates the Gaussian Blur and stores the result on the height map given
void GBlur::GaussianBlur(HFPointer inHeightMap, const int inRadius, BuildContext &inContext, const float inWeight)
    int pixels_on_row = 1+(inRadius*2);
    int height = inHeightMap->h();
    int width = inHeightMap->w();

    HFPointer temp_smap   = GetNewHF(inContext);    ///< Temporary map used for storing intermediate results

    float* gaussian_kernel = ComputeGaussianKernel(inRadius,inWeight); ///< Compute our gaussian 1d kernel

    float* pheight_map = inHeightMap->GetDataPtr();             ///< Pointer to current map
    float* ptemp_map = temp_smap->GetDataPtr();                 ///< Pointer to intermediate map
    int current = 0;                                            ///< Helps keep track of where we are

    // Do a one dimensional blur in the y direction
    // We use the temp map to find the horizontally blurred pixels
    // These are then used to compute a first blur pass and stored in the original map
    float* out = ptemp_map;
    int height_clamp = height - 1;
    int width_clamp = width - 1;

    for(int y=0;y<height; y++)
        int row = y*width;              ///< Specifies the current row

        for(int x=0;x<width; x++)
            float blurred_value = 0.0f;
            for(int xoffset = 0; xoffset < pixels_on_row; xoffset++)
                // Clamp x index
                int sx = iClamp(0,width_clamp,(x-inRadius)+xoffset);
                // Calculate newly blurred value
                blurred_value += pheight_map[row+sx]*gaussian_kernel[xoffset];

            // Set our calculated value to our temp map
            *out++ = blurred_value;

            // Increment where we are

        inContext.ReportDeviceProgress(this, y/2 ,inHeightMap->h());

    // Used for showing progress
    int half_height(inHeightMap->h() / 2);

    // Do a one dimensional blur in the x direction
    // We use the temp map to find the horizontally blurred pixels
    // These are then used to compute a second blur pass and stored in the original map
    out = pheight_map;
    for(int y=0;y<height; y++)
        for(int x=0;x<width; x++)
            float blurred_value = 0.0f;
            for(int yoffset = 0; yoffset < pixels_on_row; yoffset++)
                // Clamp our offset in y
                int sy = iClamp(0,height_clamp,(y-inRadius)+yoffset);

                // Calculate blurred value
                blurred_value += (ptemp_map[(sy * width) + x]*gaussian_kernel[yoffset]);

            // Set the original height map value to our computed value
            // This is the actual blurred value (previously scaled)
            *out++ = blurred_value;

            // increment counter and report progress to wm

        inContext.ReportDeviceProgress(this, half_height + (y/2),inHeightMap->h());

    // Release

Houdini Triangle Grid using HOM

Thought I’d post my first little (of many) experiments.
A piece of code that helped me create the header image in Houdini 11.
Sort of a hello_site example.

I plan on posting something every week or so. As in a tutorial or small experiment that helped me achieve a certain task or goal. At work or for fun.

To get this snippet of code working, create a Python geometry operator in Houdini: “file” > “new operator type”. Select: “Python Type” and “Geometry Operator”. Give the new operator a name and location on disk (default is your Houdini home / “otl”

You should be presented with a screen asking you to create new parameters and place your code.

For this piece of code to work, 4 parameters are required:

  • name: width, type:  int, default: 20
  • name: slope, type:  int, default: 1
  • name: size,   type:  float, default: 1
  • name: attr_name, type: string, default: pscale

On their own these parameters do absolutely nothing. We need to bind them in our code. See the example below on how to achieve this. All code is commented and can be placed directly in the python operator.

After pasting or writing the code, click on “Accept”. You should now be able to tab and find the operator. The “width” defines the amount of base points, the   “slope” defines the triangle angle slope, the “size” specifies the element size (or space between points). The attributre “pscale” is created and set for each point. When copying a cube on to every point, the cube should have the correct size because of the pscale attribute.

# This code is called when instances of this SOP cook.
geo = hou.pwd().geometry()

# Evaluate our parameters
width = hou.pwd().evalParm("width")
slope = hou.pwd().evalParm("slope")
size  = hou.pwd().evalParm("size")
attrib_name = hou.pwd().evalParm("attr_name")

# Initialize base numbers
height = 0
offset = 0

# Add point attribute
attrib = geo.addAttrib(hou.attribType.Point, attrib_name, 1.0)

# Compensate for slope
width_offset  = 2*slope
slope_offset  = 1*slope

# Compensate for center
center_offset = float(width) / 2.0

# Calculate Point Grid
while width > 0:
    for i in range(width):

        # Create the point
        point = geo.createPoint()

        # Get the x position
        x_position = float(offset+i) - center_offset

        # Set the point in space

        # Add point attribute
        point.setAttribValue(attrib, size)
    # Compensate the position every iteration
    width -= width_offset
    height += 1
    offset += slope_offset