What is the most powerful thing in the universe? Perhaps nobody knows, but in our world it might either be stokes theorem or the Laplace operator. It is shocking how many algorithms we end up dealing with the Laplace operator in its core steps. Chances are good that whatever you do with the knowledge of these tutorials will require the Laplace operator in some form.

In this tutorial, we shall build the cotan-Laplace operator to compute a minimal surface given a fixed boundary. Since this is a technical tutorial we will leave out the theoretical reasoning behind this. However, we will show you how to build the sparse matrix since it requires detailed technical knowledge of the relationships between points, vertices and primitives to build a sparse matrix using scipy.

##### The Input Geometry

To compute a minimal surface we need a starting surface with the topology that we shall use. In our case we will keep things simple and use a disc mesh that we will deform into some random shape. Here we display how our shape was made:

With the following code in the deformation:

1 2 3 4 5 6 7 8 |
// get distance to self float r = sqrt(@P.x*@P.x + @P.y*@P.y+ @P.z*@P.z); // and do some polar coordinate transform if(r > 0.){ float phi = atan2(@P.x,@P.z); @P.y=(sin(phi*5)*0.4+1)*log(r+1); } |

That results in this LA tostada shape:

We can also start with other shapes such as the two following ones:

# The Surface Minimization Algorithm

Now we get to the complicated part which we will guide you through.

##### Iteration Control

Our method is an iteration that can be implemented using an Block Begin – Block End pair of nodes. The complicated part of the algorithm will be put in a subnet in between that we call “Reduce Surface Area”.

Note that we can use a clever trick to quickly edit the number of iterations. In the block end node, simply replace the “Iterations” field with **$****F **so that you can edit this field using your keyboard from anyplace in Houdini just as you would normally edit the frame number.

Inside the Subnet we will have a network looking like this:

##### Python Preparations

At first we remesh the node and then make some preparations for python. We create attributes for the point Ids and for the positions.

1 2 3 4 5 6 7 8 9 10 11 12 |
// run over points: // here we store floats because it // really helps the python computations // define function g @gx = @P.x; @gy = @P.y; @gz = @P.z; // define function f @fx = @P.x; @fy = @P.y; @fz = @P.z; |

1 2 3 |
// run over points: // and enumerate pts for python i@id = @ptnum; |

Now we have to come up with a more complicate enumeration scheme. Our sparse Laplace matrix will have the number of rows and columns equal to the number of points and values will be assigned to the diagonals but also the $i$th row and $j$th column entry will be written into whenever there is an edge connecting the $i$th point with the $j$th point. How can we collect these values fast? Here is one way to do it: for every triangle collect the point ids of the 3 points around it and store them into a attribute vector resting on the primitive. Not just one vector, but two vectors as computed by the following code.

1 2 3 4 5 6 7 8 9 10 11 12 |
// record point ids for sparse matrix creation later int he = primhedge(0,@primnum); // get triangle point ids int p1 = hedge_srcpoint(0,he); int p2 = hedge_dstpoint(0,he); int p3 = hedge_presrcpoint(0,he); // store sequence for row and col // ( pt1->pt2 is one edge, etc... ) v@start = set(p1,p2,p3); v@end = set(p2,p3,p1); |

**Explanation:** The $k$th entry ( $k\in\{1,2,3\}$ ) of the **start** vector indicates from where an edge starts while the $k$th entry of the **end** vector indicates where that exact same edge ends. Now our code looks like this:

##### Marking the Boundary

Next we want to label each boundary point using an attribute. A magical Houdini trick to archive that involves the **group node**. This allows us to name a group for a specific selection of conditions (or manual selection with the mouse). In our case, we can group together the *edges* and mark the option *unshared edges*. Boundary edges are unshared edges and thus the right points are grouped together in a group that we name *boundary*.

Next we can add a point wrangle that only runs on the group: *boundary*. This way we can set an integer boundary attribute equal to 1 only for boundary points.

1 2 |
// mark the boundary with an attribute i@boundary = 1; |

Next we have to set the function g equal to zero everywhere but on the boundary. This is done in favor of the algorithm that we use and is done with the following point wrangle code:

1 2 3 4 5 |
// run over points: // set g to zero off the boundary @gx*=i@boundary; @gy*=i@boundary; @gz*=i@boundary; |

##### The Making of the Sparse Laplace Matrix

Now we will focus on creating the Laplace matrix. Sparse matrices require 4 details to be constructed.

- The dimension of the matrix: in our case the number of points.
- An array of row indices of non zero entries.
- An array of column indices of non zero entries corresponding to the above row indices.
- An array of values corresponding to the row/column index pair.

Let us try to compute the values for the matrix. The first thing we have to do is to compute edge lengths. We do this by iterating over the vertices as the number of vertices is equal to the number of half edges and Houdini does not iterate over the edges.

1 2 3 4 5 6 7 8 9 10 11 12 |
// run over every vertex: // grab the length of every half edge // get edge int halfEdge = vertexhedge(0,@vtxnum); // get position of endpoints vector P1 = attrib(0,'point','P',hedge_srcpoint(0,halfEdge)); vector P2 = attrib(0,'point','P',hedge_dstpoint(0,halfEdge)); // store distance f@edgelength = length(P1-P2); |

Next we compute the cotan values on every vertex.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 |
// run over each vertex: //compute the cotan of each vertex // make a helper function for the cotan float getCosineAlpha(float a,b,c) { float d = b * b + c * c - a * a; float e = 2 * b * c; return d / e; } // define the cotangent float getCotangentAlpha(float a,b,c) { float cosalpha = getCosineAlpha(a, b, c); float d = 1 - cosalpha * cosalpha; return cosalpha / sqrt(d); } // get edge, start vertex and end vertex int this = @vtxnum; int he = vertexhedge(0,this); int next = hedge_dstvertex(0,he); int prev = hedge_presrcvertex(0,he); // read the 3 edgelenght attachet to them float a = attrib(0,'vertex','edgelength',this); float b = attrib(0,'vertex','edgelength',next); float c = attrib(0,'vertex','edgelength',prev); // and store the cotan in the vertex. @ctan = getCotangentAlpha(a,b,c); |

Now that we have all the cotan weights ready, we only need to collect them in a reasonable way for python to grab them. For the off-diagonal entries we will grab the 3 cotan weights around each triangle and store them in one vector inside the triangle. Notice that we store the edge related cotan values inside the attribute vector in the exact same order as we have stored the *start* and *end* vector point indices in order to correctly relate the edges with their cotan value in the sparse matrix later.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 |
// run over all primitives (triangles) // compute the cotan weights of each edge in the triangle using the cotan values float getCotangentWeight(int geo; int he){ // grab "other" edge to the other face int opphe = hedge_nextequiv(geo,he); // get the two source vertices int v1 = hedge_srcvertex(geo,he); int v2 = hedge_srcvertex(geo,opphe); // read the boundary label of the points of them int bd1 = attrib(geo,'point','boundary',vertexpoint(0,v1)); int bd2 = attrib(geo,'point','boundary',vertexpoint(0,v2)); if(bd1 == 1){ return 0; // ignore boundary }else{ // the startvertex determines the row the weight is linked to. // for boundary vertices all row entries up to the diagonal // entry shall vanish float cotanleft = attrib(geo,'vertex','ctan',v1); float cotanright = attrib(geo,'vertex','ctan',v2); return -(cotanleft+cotanright)/2.; } } // grab the 3 edges and their cotan weights. // NOTE: the order is the same as in the "start","end" enumeration int he = primhedge(0,@primnum); float c1 = getCotangentWeight(0,he); he = hedge_next(0,he); float c2 = getCotangentWeight(0,he); he = hedge_next(0,he); float c3 = getCotangentWeight(0,he); // store the weights in the same order in // an attribute vector resting on the triangle v@omega = set(c1,c2,c3); |

The diagonal entries of the cotan-Laplace operator depend on all other entries in the row/column and we have one diagonal entry per point. This is why for the diagonal entries we will grab a point wrangle node and visit all the neighbours of each point to sum up the cotan weights.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 |
// run over each point: // handle diagonal entries of the laplace operator if(attrib(0,'point','boundary',@ptnum) == 1){ // diagonal entries corresponding to boundary values shall be 1 f@omega = 1.; }else{ // else sum the values on the vertex/edge star // initial edge int he0 = pointhedge(0,@ptnum); // other edges for the do while loop int he = he0; int opphe; // opposite half edge // collect the sum f@omega = 0; // cycle through every edge around a point until you get home do{ // advamce the opphe opphe = hedge_nextequiv(0,he); // collect weights of Laplacian @omega += attrib(0,'vertex','ctan',hedge_srcvertex(0,he))/2.; @omega += attrib(0,'vertex','ctan',hedge_srcvertex(0,opphe))/2.; // advance the half edge he = pointhedgenext(0,he); }while(he!=he0 & he!=-1); } |

And at last we have everything ready to really build the cotan-Laplace matrix inside python. In the next python node we will collect all the needed attributes from the points and the primitives that we have just created. For the Off-diagonal part of the sparse matrix creation we grab the indices from the *start *and *end* vectors and the corresponding values from the *omega* vector on the primitives. For the diagonal entries we grab the point indices with the attributes we just saved on the points.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 |
from scipy.sparse import csr_matrix import numpy as np import math # this code only runs once ### get geometry node = hou.pwd() geo = node.geometry() # get number of points n = len(geo.points()) ### get point and prim attributes for spase matrix creations pointid = np.array(geo.pointIntAttribValues("id")) start = np.array(geo.primIntAttribValues("start")) end = np.array(geo.primIntAttribValues("end")) # diagonal entries stored on points omegaii = np.array(geo.pointFloatAttribValues("omega")) # non diagonals stored through edges listed in primitives omegaij = np.array(geo.primFloatAttribValues("omega")) ### prepare matrix # the sparse matrix row-col references are made through the point indices # only non-zero matrix elements: # first the diagonal entries, then the entries corresponding to edges # we connect these index pairs into two long arrays rowids = np.append(pointid,start) colids = np.append(pointid,end) # values of Laplacian and mass corresponding to (row,col) specified # above directly come the omega values valueslaplacian = np.append(omegaii,omegaij) ### initialize matrices laplace = csr_matrix((valueslaplacian,(rowids,colids)),(n,n)) # cache matrix locally for the next python node node.setCachedUserData("laplace",laplace) |

##### Solving the Poisson Equation

Next we will show how the Poisson equation $\Delta f=g$ can be solved inside a python node.

First we will read the function g into python (we had stored it as float attributes *gx,gy,gz* and we then cache them into python.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
import numpy # read this node node = hou.pwd() # read the function g values gx = numpy.array(node.geometry().pointFloatAttribValues("gx")) gy = numpy.array(node.geometry().pointFloatAttribValues("gy")) gz = numpy.array(node.geometry().pointFloatAttribValues("gz")) # here we cache the function g locally for the next node to access node.setCachedUserData("gx", gx) node.setCachedUserData("gy", gy) node.setCachedUserData("gz", gz) |

And next we take another python node, read all the inputs g and the Laplace matrix in and solve the Poisson equation.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 |
from scipy.sparse.linalg import spsolve import numpy # grab the node node = hou.pwd() geo = node.geometry() # extract the laplace matrix laplace = node.inputs()[0].cachedUserData("laplace") # extract the function g gx = node.inputs()[1].cachedUserData("gx") gy = node.inputs()[1].cachedUserData("gy") gz = node.inputs()[1].cachedUserData("gz") # solve the poisson problem # use spsolve since our matric is symmetric and pos. def. fx = spsolve(laplace,gx) fy = spsolve(laplace,gy) fz = spsolve(laplace,gz) # store solution into attribute for VEX code geo.setPointFloatAttribValuesFromString("fx",fx.astype(numpy.float32)) geo.setPointFloatAttribValuesFromString("fy",fy.astype(numpy.float32)) geo.setPointFloatAttribValuesFromString("fz",fz.astype(numpy.float32)) |

And at very last we copy the computed values into the point positions using a point wrangle node.

1 2 3 4 |
// copy the attribute values back to the positions @P.x=@fx; @P.y=@fy; @P.z=@fz; |

# Displaying the Results

We want to display the results using soap films with highlighted boundary. Here we present a simple way how to make that possible.

To highlight the boundary we label every boundary point again and then use that information in a detail attribute wrangle node in order to loop over each vertex (thus over each edge) and add a line segment to that edge if it is a boundary edge. The fuse node afterwards then takes care of fusing two points from different line segments into one point. This way we actually have a complete polygonal line for every boundary component. Here is the code for the creation of the boundary lines. For the textures please look at the texture tutorial.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 |
// this code runs only once (detail) // here we run through every vertex to run over each edge. // if an edge is between two boundary points we create // a polyline between them // the redundant points will be merged by the fuse node // after this node. //in: geometry int, 2 point ids //out: a new line in the geometry int addline(const int geo, p1, p2){ int prim = addprim(geo, "polyline"); addvertex(geo,prim,p1); addvertex(geo,prim,p2); return prim; } // determine references int srcGeo = 1; int dstGeo = 0; // run over every vertex int numvtx = nvertices(srcGeo); for(int vtxnum=0;vtxnum<numvtx;vtxnum++){ // grab edge int edge = vertexhedge(srcGeo,vtxnum); // grab points int start = hedge_srcpoint(srcGeo,edge); int end = hedge_dstpoint(srcGeo,edge); // check if both at boundary int startIsBd = attrib(srcGeo,"point","boundary",start); int endIsBd = attrib(srcGeo,"point","boundary",end); // create polyline there if both are boundary if( 1 && startIsBd==1 && endIsBd ==1){ // add two points first vector Pstart = attrib(srcGeo,"point","P",start); vector Pend = attrib(srcGeo,"point","P",end); int p1 = addpoint(dstGeo,Pstart); int p2 = addpoint(dstGeo,Pend); // then add the line addline(dstGeo,p1,p2); } } |

The results looks as follows:

Note that if you make the tube too long the minimal surface has to split the mesh into two discs. However, our implementation does not yet handle such geometry splittings.

Just to clear things up, this is how our final network looks like: