It seems like I’m about to head into the dark territory of F# “quotations”, which apparently let you play compiler without all that tedious parsing. (Sorry for the absurdly shallow interpretation, but I haven’t played with them much yet.) Anyway, before I get into that, I thought I’d write up my progress on getting nicely rendered space curves.

### Segmented Cylinders

Last post, I’d gotten a unit cylinder to render halfway decently, after some compromises (and throwing more triangles at the thing). After that, I had to translate and rotate them where I wanted them. That should have been a trivial step (“I have a unit vector along the *x* axis, please map it over here”), but turned out to be more trouble than it should have been. (I won’t tell the full story, but here’s a quick runthrough.) The Matrix.CreateWorld method looked very promising indeed, but kept squashing out any *z* component (for reasons which became clear later; at the time I was passing the z unit vector as my “up vector”). After some hunting, and reading far more than I ever wanted to about quaternions, I fell back to a simpler approach: I want to overlay the 0-1 line on the *x* axis onto some line *p0* to *p1*. Clearly, I need to translate by *p0* and scale by |*p1-p0*| (that’s supposed to be the length of the line, if my pseudo-typography isn’t coming through). In almost all cases, I also need to rotate; from basic vector math, the axis of rotation will be the cross product *i* × (*p1-p0*) (I can’t get the little hat over the i, but that’s the unit vector in the *x* direction). The amount of the rotation will be the angle between the vectors: *acos(i • (p1-p0))*. (I’m particularly fond of UK-based sites, since they tend to refer to “maths”. Juvenile, I know.) This all works fine except for the case when *p0* and *p1* actually lie on the x axis already; in that case, the cross product goes to all zeros, which was what was flattening out my matrices from CreateWorld (at least I think so; I never went back to test the theory). So I ended up with a scale, a conditional rotation, and then a translation. I was sad about the conditional — even though modern shaders can handle conditionals, it still feels like a mark of shame. There’s probably some way to be clever about this; after all, **any** axis of rotation will do just fine if you’re rotating zero degrees. But I was in a hurry to get something rendering again, so I gave up for the moment.

Anyway, this was enough to get me from oddly-lit unit cylinders all the way up to this:

In this case, I’m rendering the first Hermite basis function, since I had code lying around to calculate it. With 20 samples, the curve looks quite smooth (all things considered), but it’s discontinuous, and gets into worse trouble near the ends. (In areas of high curvature, the abruptly changing normals at lines of overlap look distinctly articulated, rather than like a smooth curve.) So the next step was to get at least C0 continuity, and try to do something about the normals.

### Beveled Cylinders

At this point I was on a roll with matrices. Nothing feels as much like a triumph as finally succeeding at a seemingly trivial task that took you all Thanksgiving afternoon (between moments of wincing at the Lions). I decided I wanted the cylinders to join along the plane which passes through the actual endpoint *p_n*, and whose normal is the average of the incoming and outgoing vectors, as seen in the following (atrocious) Paint diagram:

The plan is to virtually slice the cylinders along that plane. (Note that the average of the incoming and outgoing position vectors was chosen so that both cylinders would be sliced at the same angle; otherwise the diagonal sections wouldn’t match up.) This means that approximately half the vertices will be placed farther along the position vector than they were in the sandworm version, and half will fall short of their positions there. This approach preserves the cross-section of the cylinder — at any point, a cut normal to the long axis of the cylinder would show a perfect circle or a segment of one — but potentially allows arbitrary position error at the endpoints. In particular, if two segments meet at a very narrow angle, the outermost points may have to be projected a very long way in space to get to that common plane. (There are alternative strategies, as I consider later on; first I just wanted to get **one** of them to work. Besides, meeting at such a narrow angle is a pretty good indicator of insufficient tessellation in the first place; either it is just a nasty curve, or we’re sampling too sparsely in an area of high curvature.)

I should stop at this point to note the obvious: I could quite simply guarantee C(0) continuity by emitting multiple cylinder segments as a continuous strip, rather than transforming a single cylinder to a series of locations. I didn’t want to start with that approach, however, because I wanted to have a good justification for why both segments should share those endpoints, rather than simply force majeure. I also liked being able to play with different matrices for each end of the cylinder, which would have been trickier if they were all passed in a single batch. (Not insurmountably so, but I wanted to focus on the geometry, not on instancing techniques.)

Anyway, I had a clear description of where I wanted the cylinder endpoints to end up. Now the trick was to get them there. The two ends are necessarily independent, since the meeting angles at opposite ends of a segment have nothing naturally in common (although now that I think about it, that would be an interesting tessellation heuristic; even then it would probably be an upper limit, rather than a fixed value). I’d just been reading about geometry instancing, trying to figure out how I was going to deal with tens of thousands of these cylinders running around, so the natural approach was to pass an array or arrays of shader parameters. The x coordinates of the cylinder vertices are all 0 or 1, so I can use them to index the arrays, and then be freely able to specify as many parameters as I need for each end independently. (I was a little surprised at using a floating point coordinate as an array index, but it works fine, and seems to be the standard approach.)

So now the ends of the cylinders are decoupled. I’m trying to project the vertices of each end onto a (different) plane, so the natural thing to do is to pass an array of two matrices, one for each end. Matrices are good at projections, after all. Trouble is, the projection that comes to mind would all project the points along the normal to the plane of projection — an orthographic projection operates this way. But I want to project *along the cylinder* and intersect the plane wherever that happens, rather than taking the shortest path. (Note that projecting along the plane’s normal would end up distorting the cross section of the former “cylinder”; in most cases I think it would turn into an extruded ellipse instead of a circle.) But just as I started trying to work out the math from scratch (and old linear algebra texts), something made me look at the static constructors available on the Matrix class again. Lo and behold, there’s the CreateShadow method: “Creates a Matrix that flattens geometry into a specified Plane as if casting a shadow from a specified light source.” Using the long axis of the cylinder as my “light source” direction, and the plane described above as the projection plane, this fits the bill perfectly! (There are some follow-on problems, but first pretty pictures from when I was still happy this was working.)

This looks pretty good for only 10 samples. In a larger version you can see the specular highlight “bleed” a little along the knuckle of the vertex in the bottom right, but this is at least partly due to the abrupt change in the angle per pixel rate, as I talked about in my previous post.

Here is a wireframe version, where you can quite clearly see the endcaps of the cylinders being extended and contracted around the sample points:

All is not altogether wine and roses, however. There are fundamental problems when the local radius of the curve is tighter than the radius of the cylinder: the “tube” starts to overlap itself and develop strange buckles and creases. There also seems to be an additional problem, because I’m not sure that that by itself would explain all the oddities in this badly-formed sine wave:

Sad to say, this test also demonstrates the problems of arbitrary positional error which I mentioned earlier; the massive “hoods” at either end don’t do a great job of reflecting the shape of the curve (and have strangely messed-up normals). I think part of the problem here is that I locked the free ends of each curve (where there’s only one segment, hence no average plane) to the vertical axis. This looks very pretty on curves like the Hermite function above, and makes for good screenshots of the [0..1] region, but I suspect from this picture that I’d do better to use the one segment I do have, rather than an arbitrary fixed value.

### An Afterword: Normals

An issue with this approach is how it handles normals. I intentionally used the average of the as-tessellated position vectors, *not* the analytical tangent, to choose the plane which is the locus of the endpoints. The goal was to improve quality when the as-sampled and analytical tangents diverge. First, consider a sine wave which is only sampled at the ±1 points. In each case, the analytical tangent is zero. However, we necessarily have to travel vertically to get to the next *sampled* point. Using the analytical tangent would tend to produce a shearing effect on the cylinder, while the as-sampled approach gives a nicely angled corner. As a second example, consider the same sine wave, but only sampled when it crosses the axis. Every sample is zero, but using the analytical tangent would have the seams slanting back and forth to no particular gain.

The trouble, of course, is that my normals are getting flattened into the same plane as the vertices; this effectively gives me OpenGL-style “computed normals” as the average of the adjoining segments. (Note that I’m not 100% sure *where* my normals are actually ending up; this description agrees with what they look like, but I’ve been burned by that before. The point of confusion is what effect I’m producing by transforming my (non-homogenous) normal vectors by the upper-left 3×3 matrix, after multiplying in the shadow matrix. It sure looks like they’re getting projected into the plane, but I don’t want to swear to it.)

To go a step further, I would want to calculate how much the analytical tangent at the point deviates from the as-sampled tangent, and use that correction factor as the basis for an additional rotation of the normals (potentially rotating back out of the shared plane). I have not yet pursued this approach because the whole strategy is assuming the availability of the previous and next point in the tessellation, and I haven’t yet decided how well that will scale to serious instancing. (I did play with it, passing the un-shadow-matrixed translate/scale/rotate matrix as a third matrix in the array, and got some tidy-looking curves looking like they were made from segments of PVC pipe, but that wasn’t really the goal.)

However, before I get into trying to solve that puzzle, I want better control over what I’m doing. I should probably shell out the $99 for nVidia’s shader debugger, but I’m much more interested in being able to write my shader code in F#, to evaluate on either the GPU or the host processor. I suspect that this exercise will be a good chance to get back to my functional programming roots, give me a pretty deep understanding of F# mechanics, and also serve as a good introduction to the nuances of shader evaluation.

Having some spare time on my hands right now, I came to read your bit on functional programming – which I heart by the way – but I had to comment on this entry to agree that “maths” is really the most awesome usage and should be applied whenever possible.

Comment by Vanessa — March 8, 2010 @ 3:20 pm