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

'''
Globals______________________________________________________________________________________________________
'''

# 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

'''
Attributes___________________________________________________________________________________________________
'''

# Add the number of connections attribute
num_connections_attrib = geo.findPointAttrib(num_connections_name)
if num_connections_attrib is None:

# Add the number of associations attribute
num_associations_attrib = geo.findPointAttrib(num_associations_name)
if num_associations_attrib is None:

# Add the point resolved attribute
resolved_attrib = geo.findPointAttrib(unresolved_name)
if resolved_attrib is None:

resolved_prim_attrib = geo.findPrimAttrib(unresolved_name)
if resolved_prim_attrib is None:

'''
Methods______________________________________________________________________________________________________
'''

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

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

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 group.name() == inGroupName:
if len(group.points()) > 2:
pgroup = []
for point in group.points():                            # Search for the point in the original points and return those!
pgroup.append(points[point.number()])
return pgroup
else:
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 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):
"""
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])
else:
edge = Edge(inPoints[pnum], inPoints)

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 group.name():
# 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, group_points)
active_edges.append(edge)
boundary_edges.append(edge)

# 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:
rlist.append(attr)
else:
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)

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))
active_edges.remove(edge)
continue

if edge in unresolved_edges:                                   #if in unresolved edges, remove
Console("edge was in unresolved edges: %s" % str(edge))
unresolved_edges.remove(edge)
continue

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))
vector.normalized()
prim_point_normals.append(vector)

# combine normals
accum = hou.Vector3(0,0,0)
for n in prim_point_normals:
for i in range(3):
accum[i]+=n[i]

# 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()
if posx < xmin:
xmin = posx
order_array = point

inPoints.remove(order_array)

for point in inPoints:                  # find second point
posz = point.position()
if posz < zmin:
zmin = posz
order_array = point

inPoints.remove(order_array)         # append third point
order_array = inPoints

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()
poly.setIsClosed(True)

if inEdge.pointOne.isEdgePoint and inEdge.pointTwo.isEdgePoint:
order_array = [inEdge.pointOne, inEdge.pointTwo, inPoint]
else:
order_array = SortVertices([inEdge.pointOne, inEdge.pointTwo, inPoint])

# add the poly as a TriPoly to our list of polygons
bound_polygons.append(TriPoly(order_array))

# 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
UpdateActiveEdges(edge_list)

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

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

p2_dist = inDictTwo[key_point]                                                   # get the point two distance
p1_dist = inDictOne[key_point]                                                   # get the point one distance

if point_prior == "minimal":                                                        # choose sampling method
dist = abs(p1_dist - p2_dist)
else:
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

"""
Main------------------------------------------------------------------------------------------
"""

# 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:
CreateSeedEdges(GetSeedPoints(edge_group_name))
else:
CreatePredefinedSeedEdges(edge_group_name)

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

current_edge = active_edges
closest_point = GetEdgeNeighbour(current_edge)
if closest_point != None:
CreatePolygon(current_edge, closest_point)
else:
unresolved_edges.append(current_edge)
active_edges.remove(current_edge)
count += 1

Console("\n")

Finalize(unresolved_group)

# 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()
edge_point.setPosition(edge_point_position)

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

attrib = geo.addAttrib(hou.attribType.Point, "edge_normal", (0.0,0.0,0.0), transform_as_normal=True)

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

p1.setPosition(pp1)
p2.setPosition(pp2)
p3.setPosition(pp3)

point_list = [p1,p2,p3]

# Create Tri poly
poly = geo.createPolygon()
poly.setIsClosed(True)

# 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]
else:
point_one = point_list[-1:]
point_two = point_list

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)
{
float* gaussian_kernel = (float*)malloc(mem_amount*sizeof(float));

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

sum+=v;
r++;
}

// 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 height = inHeightMap->h();
int width = inHeightMap->w();

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

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

// 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
current++;
}

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

// 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
current++;
inContext.ReportDeviceProgress(this,current/2,inHeightMap->area());
}

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

// Release
free(gaussian_kernel);
gaussian_kernel=NULL;
temp_smap->Clear();
temp_smap->releaseAccessible();
}

# Houdini Grid Based Fluids And Dynamics

A while back I got an invite to do a two day seminar in Seoul.

Teaching at the university of Seoul about Houdini’s Dynamics environment and recursive functionality. After waiting for over a year I find it’s about time to release the training material.

Below you will find a link to a rar package containing a pdf, example files and demo otls.

The material consists of an introduction in to recursiveness in Houdini, Houdini’s Dynamics environment, Houdini Volume primitives and fluid solvers. Starting simple and ending with building a custom fluid solver incorporating a simple Dynamic Resize model and simple volume shape operators.

Hopefully enough to get anyone going within Houdini’s dynamic environment. But most importantly, to create awareness regarding the full set of possibilities this environment has to offer. Never again should you wonder what this magic shelf button does.

Also note that some of the .hip files might have incomplete geometry. Creating a simple new geo reference should suffice. I will try to get the files up to date as soon as possible.   # 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

# 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
point.setPosition((x_position*size,height*size,0))