Showing posts with label Space Colonization. Show all posts
Showing posts with label Space Colonization. Show all posts

Saturday, 10 January 2015

Space Colonization Algorithm part 3

Holidays are over unfortunately which means time is more constraint. Most of what is in this part I did on the last days of my holidays and the changes to the source code have already been added to the GitHub repository.

First I had to make a few improvements to the tree optimisation algorithm. I had actually made a few typos in some of the vector classes and it wasn't calculating things correctly. Once that was all ready I could tweak the optimisation algorithm slightly and the result is great. The only noteworthy tweak is that I calculate and keep the direction vector as long as I'm merging the same section. This results in very slow bends still being optimised much better.

OpenGL 3/4

The most work went into changing over to using OpenGL 3/4. As I'm using a tessellation shader step for automatic LOD (more on this later on) the program now requires OpenGL 4 capable hardware. If that isn't what you have you could remove the tessellation logic and enhance the mesh generator to add more detail.

I'm not going into much detail about this step, I might do a separate (set of) blog post(s) on setting up OpenGL 3 as a render environment. But I'll go through the basics.

As I'm developing on a Mac most of the fixed rendering pipeline is gone. It basically means you're doing everything in shaders yourself. I've had to add a matrix class (only basic functions supported for now) and add a few things to the vector classes to make things easier. Again I'm not switching to an existing maths library as it is easy to just build these classes with just what you need but you could get a head start by using something like GLM.

I did add the excellent stb_image.h library into the mix to load textures.
Also note that in constructing the GLFW window a number of hints have been added that are required to enable OpenGL 3+ (and get you to the point of no return):

// make sure we're using OpenGL 3
glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 3);
glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 2);
glfwWindowHint(GLFW_OPENGL_FORWARD_COMPAT, GL_TRUE);
glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE);

Shaders

First we need to implement shaders now that we no longer have our fixed rendering pipeline. The shaders are all inline which makes them a little harder to read but its just easier to deploy them that way. If I had a larger project I would create them as files on disk and just load them in (in fact the game engine I'm working on comes with a complete precompiler for shaders).
There are two shaders:
- "simpleshader" which basically passes through a color
- "treeshader" which renders our final mesh and has a tessellation shader step

Lets have a quick look at my simple shader as it explains a few things on how this stuff works and is used to render what we used to do in OpenGL 1.

-- vertex shader
#version 330

uniform mat4 mvp;
layout (location=0) in vec3 vertices;

void main() {
  vec4 V = vec4(vertices, 1.0);
  gl_Position = mvp * V;
}
-- fragment shader
#version 330

uniform vec4 color;
out vec4 fragcolor;

void main() {
  fragcolor = color;
}

Both shaders start with "#version 330" which simply tells OpenGL what feature set we want. Our tree shader is set to "410 core" as our tessellation shader isn't supported else wise. After 330 basically new stuff has been added but before there are difference in dialect. For instance we used to have build in variables for the color and our model/view/projection matrix but no more.

Our vertex shader is where most of our logic resides. Let's have a look at each line of code a little closer starting with our 3rd line:
uniform mat4 mvp;
The first word on this line, uniform, tells the shader this is a variable that will be set from outside of the shader but is otherwise immutable. Our variable in this case is a 4x4 matrix called mvp (model/view/projection matrix).
In our rendering code we see that we calculate our model/view/projection matrix and then give it to our shader using the setMat4Uniform method on our shader class which in turn calls glUniformMatrix4fv which is our OpenGL command for copying the matrix.
layout (location=0) in vec3 vertices;
This one is a bit more tricky. "layout (location=0)" tells OpenGL this is our "0" attribute in our vertex buffer object. We bind this by calling glVertexAttribPointer. Before OpenGL 3 this already existed but was somewhat hidden and hardcoded. Also before OGL3 you had to call this every time before you rendered your mesh while it is now stored in our vertex array object (more on this later).
To make a long story short, this line of code makes our vertices available to us in our shader and automatically refers to the correct vertex being handled in our shader.
void main() {
Indicates the start of our shader program
vec4 V = vec4(vertices, 1.0);
Simply turns our 3D vertice into a 4D vector
gl_Position = mvp * V;
Applies our matrix to our vertex so that we end up with the correct coordinates on screen and in our depth buffer.

Our fragment shader is much simpler
uniform vec4 color;
Provides us with a way to set the color we want to render with from outside of our shader.
out vec4 fragcolor;
Defines that the output of our shader is a RGBA color variabled called fragcolor.
fragcolor = color
Simply assigns our input color to our output.

The shaders are compiled once and then simply used by calling glUseProgram.
For compiling the shaders have a look at the shader class in the source code.

Vertex Array Object 

The second change we must do is that we need to bind all our state into a Vertex Array Object or VAO for short. A VAO is little more then a container of state, before VAOs were added to OpenGL you had to first bind all the buffers and related state before rendering the object.
VAOs simply added some convenience allowing your to do all the binding calls once and then simply make the VAO active.
In OpenGL 3 they've made the use of a VAO mandatory which can be a bit of a nuisance.

If you're doing the right thing before you render anything, you would create a VAO for your 3D object, bind the buffers, load them with data, and then reuse the VAO when you need it during rendering.
In my tree application I do some of this setup in my render loop and sometimes even repeat it needlessly. This is purely due to wanting to keep showing each step of the algorithm and I wouldn't recommend it as the right approach:)

Lastly, the VAO binds Vertex Buffer Objects to itself and when activated reloads that state. The contents of those buffers is another matter entirely and you don't need to have a VAO active to update or set the contents of the buffers. Equally so, a buffer object is not restricted to a single VAO. If you have a mesh that is rendered using different shaders to simulate different materials you might create a VAO for each of those materials, bind the same vertex buffer that contains the vertices of your mesh but bind a different buffer that contains the indices of the faces being rendered.

Vertex Buffer Objects

Vertex Buffer Objects, or VBO for short, have been around for awhile now. They are very simply put buffers of data, generally speaking either a buffer containing the vertices of your mesh (and their normals, texture coordinates and other such things) or indices to form the faces of your mesh (triangles mostly).

As mentioned before, generally you would load up your mesh, create and load date into your VBOs, bind them to a VAO and  you'd be done, the VAO would be ready for rendering the mesh in your render loop.

For practicality I'm doing most of this logic inside of the render loop for the tree application. But thats generally not what you would do in a normal application.

There is more in relation to the OpenGL 3 conversion of the code but these 3 are the major highlights.

Mesh generation

With all that ground work done we finally come back to our tree generation algorithm. We've basically got a working model of the Space Colonization Algorithm already but now we need to turn it into a 3D mesh that can be rendered.

I'm taking a shortcut here because I've moved some of the complexity into the shader logic. This allows me to simply box up the tree and create a very simple mesh. It also means my mesh is made up of quads, not triangles.

For each vertex of my tree node I need to create 4 vertices within the correct plane. That plane is defined by taking the direction of the branch to that vertex and the direction to the next vertex and averaging the two, and using that as a normal for the plane.

Where a tree branches I create multiple sets depending on the direction I branch into.

After that I simply create a box using the 4 vertices of each of the 2 points of a node.

Creating our 4 vertices for our cross section was the tricky part and I'm only 90% happy with the solution I came up with. I simply take the cross product of my plane normal vector and an arbitrary constant vector to get a vector that is perpendicular to my normal vector all pointed roughly in the same direction. Then I rotate that vector in 90 degree steps to create the other 3 vertices. Finally I scale the vector depending on how far down the tree I am so the tree is thin at the end of its branches and thick at its roots.

It works really well unless my tree is growing in roughly the same direction as my arbitrarily chosen vector.
Also it didn't work very well with my roots so I've disabled them for now.

The last step is done within the tessellation control shader and the tessellation evaluation shader.

The control shader is basically a straight forward control shader were I take the size of each edge of my quad as its projected on screen and multiply it by a constant. Remembering that our screen coordinates are currently in the -1.0 to 1.0 range, not at screen resolution, we thus end up with a number that neatly divides our mesh into roughly equally sized triangles. The closer our tree, the more triangles.

If we take this division as is we would have a nicely highly detailed tree which would still be square. Our evaluation shader therefor applies a smoothing algorithm and we get a really nicely rounded mesh. For now I've used a technique called Phong Tessellation which works very well. It's described here much better then I ever could: http://liris.cnrs.fr/Documents/Liris-6161.pdf

Still I'm pretty happy with the end result:


Next steps

The most obvious omission is that we don't have any leaves. That will have my next focus.
The other thing that I've been planning to add is to add interface controls for creating the tree, maybe give the user the ability to "paint" the point cloud and define our starting tree shape.
What is also needed is a way to load/save our tree as we're working on it and to export it into a mesh format that can be used in other applications.

Saturday, 3 January 2015

Tree update

The last few days have kept me busy with the kids, playing to much on my XBox, trying to finish A Dance With Dragons and then some, but I did manage to put more time into my tree thingy.

I haven't had time to write up part 3 of the series and when I do it will be a long one as I've got an initial version of the code working that extrudes the tree nodes to a mesh, ported everything I had written over to using Vertex Array Buffers, Vertex Buffer Objects, programmable shaders including using tessellation shaders, etc. That last one requires Open GL 4.0 support but hey...

As my holidays end in a little over 24h it may be some time before I can do the write-up but for the time being, just a few progress screenshots:



I did upload all the changes to my github page.

Monday, 29 December 2014

Space Colonization Algorithm Part 2

One issue that I ran into a few times after finishing part one is that sometimes a branch would end with equal distance between two attraction points. When adding a new node I call the method called growBranch that adds the new vertex and node to the tree.
This logic has been enhanced in 3 small ways:
- I look up the "parent" node to make further processing easier
- I update a counter for the number of child nodes, this will allow us to vary the thickness of the tree later on.
- I calculate the dot product of the vector of the parent node and the new node to see if it is backtracking on itself (which is what happens in the aforementioned situation) and branch out in a separate direction using the cross product of the two vectors.

Calculating the dot product of two normalised vectors gives us the cosine for the angle between the two vectors, a very handy and oft used little calculation. A value of 1.0 means the vectors are in the same direction, a value of 0.0 means they are at right angles and a value of -1.0 means they are in opposite direction (which is what we're checking for).

The cross product of two vectors gives us a vector that is at right angle to the two vectors. Very handy.


Interface changes

Another small change I added is that I no longer rotate the model automatically but rotated it based on scrolling horizontally and going closer/further away using vertical scrolling. This works really well on a MacBook pro as these actions are controlled by using two fingers on the touchpad but may be unwanted using traditional mice, you may want to stray away from this.

I'm also drawing larger dots for the vertices of the tree to better show the output of the algorithm. 

Optimising the nodes

The real enhancement and the subject of this blog post is the added optimisation step.
After generating our tree (and zooming in a bit), we can see how many segments we've created:

This level of detail is overkill for what we're doing. Again our trusty dot product comes to the rescue. If we have two nodes whose vectors are pointing in (roughly) the same direction we could merge the nodes to a single node between the two end points and remove the point in the middle. We can only do this if a node has only one child node, if it has multiple branches the detail remains. We thus need to test for that situation.

void treelogic::optimiseNodes() {
  unsigned int node = 0;
  std::vector<unsigned int> children;
  
  while (node < mNodes.size()) {
    unsigned int mergeWith = 0;
    
    // we need to find out how many children we have, we can only optimise if just one is found
    children.clear();
    for (unsigned int n = node+1; n < mNodes.size(); n++) {
      if (mNodes[n].parent == node) {
        children.push_back(n);
      };
    };
    
    // only one child? check if we need to merge
    if (children.size() == 1) {
      vec3  parentVector = mVertices[mNodes[node].b] - mVertices[mNodes[node].a];
      parentVector = parentVector.normalized();      

      vec3  childVector = mVertices[mNodes[children[0]].b] - mVertices[mNodes[children[0]].a];
      childVector = childVector.normalized();
      
      // use dot product, this gives our cosine, the closer to 1.0 the more the vectors match direction
      if ((parentVector % childVector) > 0.9) {
        mergeWith = children[0];
      };
    };
    
    // and merge
    if (mergeWith!=0) {
      unsigned int eraseVertice = mNodes[node].b; // should be same as mNodes[mergeWith].a, this we'll erase..
      
      // merge our nodes...
      mNodes[mergeWith].a = mNodes[node].a;
      mNodes[mergeWith].childcount = mNodes[node].childcount;
      mNodes[mergeWith].parent = mNodes[node].parent;
      mNodes.erase(mNodes.begin() + node);
      
      // erase our vertice
      mVertices.erase(mVertices.begin() + eraseVertice);
      
      // adjust our nodes
      for (unsigned int n = 0; n < mNodes.size(); n++) {
        if (mNodes[n].parent > node)mNodes[n].parent--;
        if (mNodes[n].a > eraseVertice) mNodes[n].a--;
        if (mNodes[n].b > eraseVertice) mNodes[n].b--;
      };      
    } else {
      node++;
    };
  };
};

This removes most of our in between nodes:


At this point we can see that we've removed a little to much detail because some of our curves have gone. This is because the small angles between nodes when there are a lot of segments means we end up removing too much.  I'm still playing around to see if this can be improved.

Next up however, taking our nodes and turning it into a mesh that we can then further process. 

Space Colonization Algorithm Part 1

Hitting a roadblock with the model I was using to test my engine I found myself working on a few related projects.

The first was implementing a texture blending technique that was described here: http://www.nongnu.org/techne/research/blending/blending.pdf , I'll write some more about that at a later stage but it was fun to do and works amazingly well.

The second was that trees and vegetation were going to be needed at some point in time. Now I could simply import models of trees and other things but I wanted to find out if it would be possible to generate these and do this in a way that would allow me to create multi LOD meshes for the same tree. Searching the web you find a lot of neat tools and small tutorials.

Most seem to favour a technique called L-System fractals. Google on that term if you wish but it is basically a way to notationally describe a fractal system that can be used to generate a tree (and many other things). However there is a lot of trial and error involved and the more complex and fully grown a tree becomes, the more roadblocks your seem to hit.

Eventually, buried deep down in the search results someone mentioned a technique called "Space Colonization" which I promptly ignored due to the name suggesting it had more to do with space combat games then procedurally generating vegetation. But after it came back a few times most notably on this amazingly cool blog I starting looking into it. Eventually it lead me to the the very well written paper called "Modeling Trees with a Space Colonization Algorithm" written by Adam Runions, Brendan Lane, and Przemyslaw Prusinkiewicz.

I decided to try my hand at implementing it and see how far I'd get. It turned out to be very easy to do the basic generation of the tree. As I've not found many people showing how they've implemented the algorithm I decided to dedicate a few blog posts about it.
This first one is just about the node generation.
The next step will be to optimise the output merging branches, then we'll need to generate a model that will allow us to render the tree and finally I need to add in leaves and such.
Finally I need to add more algorithms for creating differently shaped point clouds or maybe even allow for drawing out the shape.
Not all necessarily in that order:)

Source code, Github and framework

I've uploaded and will be maintaining the project to my GitHub page. I'm working on a Macbook Pro at the moment so everything will only be tested on that platform for now. Obviously you'll need XCode but as I'm using makefiles (I'm a little old fashioned) you'll also need to install the command line tools and possible make (I can't remember, I set it all up years ago). Or off course you start a new project in your favourite IDE and just draw the source files in.

As I'm using OpenGL and GLUT is pretty useless once you go beyond OpenGL 2 (on a Mac at least) I've started using GLFW and GLEW. The GitHub repository contains the header files and static universal libraries compiled for Mac OS X. I just added those to make it easier to duplicate the repository and get started but you'd do well to get the original source code and documentation for both. I won't waste any words here on setting up and compiling these (I might be tempted to do so at some point in time) as the documentation is pretty good.

Right now I've not added a math library but just added those classes and those methods I've needed so far. This may change as the project progresses.

OpenGL 1.0

At this stage the project is using OpenGL 1.0 and worse even, using glBegin/glEnd commands. This is simply done because I needed to test the algorithm and it is the quickest way to just output stuff to screen. It won't stay so for much longer though. I haven't decided yet but it is likely I'll use tessellation shaders only available in OpenGL 4 to generate the detail for the mesh. Also as my main engine runs OpenGL 3 code I'll eventually end up changing the code to using proper vertex arrays and buffers.

Attraction Point Cloud

This is the heart of this algorithm. The basic idea behind the algorithm is that you have a cloud of points to which the tree grows. The shape of the point cloud very much dictates the overall look of the tree but inherent in the algorithm is the branching you would expect with vegetation.

The code I currently use to test the algorithm isn't very imaginative in the cloud it creates, I simply randomly generate positions of points within a hemisphere and stretch that to roughly represent the shape of a tree. How its stretched and how the points are organised can be set with parameters to the method that generates the point cloud. The samples in the aforementioned paper use much better algorithms producing better shapes and I plan to enhance this eventually and take it even further then the paper suggests.

Note that for this reason also the sample code animates the algorithm, this algorithm works very well visually to see how it works and see how changes to the parameters effect the outcome. When you use it for real obviously you'd just run the algorithm and work with the finished mesh.

The code to generate the attraction points is very simple:
void treelogic::generateAttractionPoints(unsigned int pNumOfPoints, float pOuterRadius, float pInnerRadius, float pAspect, float pOffsetY, bool pClear) {
  // Seed our randomiser
  srand (time(NULL));

  if (pClear) {
    // Clear any existing points (shouldn't be any..)
    mAttractionPoints.clear();
  };

  // Add random attraction points until we reached our goal
  for (unsigned int i = 0; i < pNumOfPoints; i++) {
    vec3 point;

    // random normalized vector for half a hemisphere
    point.x = randf();
    point.y = randf(0.0f, 1.0f);
    point.z = randf();
    point = point.normalize();

    // Scale it up to a random radius and stretch if needed
    point *= ((pOuterRadius - pInnerRadius) * randf(0.0f, 1.0f)) + pInnerRadius;
    point.y *= pAspect;
    point.y += pOffsetY;

    // and add it to our buffer
    mAttractionPoints.push_back(attractionPoint(point));
  };
};

You can call it a few times consecutively to create slightly more interesting clouds.
And it gives a point cloud like so:
All the green dots are the points that were generated.

Growing the tree

The magic happens in the iterative logic that grows the tree. Notice in the previous screenshot the brown line, that is the start of our tree. The idea is that you can start off with a basic shape to further control the end result. All the algorithm needs is at least one starting vertex at 0,0,0 which is taken care of in the constructor of the class.

For every iteration we check each point in our attraction point cloud and check which vertex within our tree lies closest to that point. If the distance is within a certain threshold (dk in the paper, pMaxDistance in the code below) we take it into account, else we ignore the attraction point for now. In this manner we find all attraction points that are closest to each vertex of the tree (these are our sets). For each vertex of the tree we calculate the average normalised vectors to these points and then branch out in that direction by a specified distance (D in the paper, pBranchSize in the code beow).
Finally, any attraction point that is closer then a specified distance (di in the paper, pCutOffDistance in the code below) to its nearest vertex is discarded.

We repeat this process until all attraction points have been discarded.

The code for our iteration looks like this:
bool treelogic::doIteration(float pMaxDistance, float pBranchSize, float pCutOffDistance, vec3 pBias) {
  unsigned int numVerts = mVertices.size(); // need to know the number of vertices at the start of our process
  unsigned int i, v;
  std::vector<float> numOfAPoints;
  std::vector<vec3> directions;
  
  // init our temporary buffers
  for (v = 0; v < numVerts; v++) {
    numOfAPoints.push_back(0.0);
    directions.push_back(vec3(0.0f, 0.0f, 0.0f));
  };
  
  // find out what our closest vertice to each attraction points is:
  i = 0;
  while (i < mAttractionPoints.size()) { // use a while loop, we'll be removing points..
    attractionPoint point = mAttractionPoints[i];
    
    // start with our current distance for our attraction point
    vec3 delta = mVertices[point.closestVertice] - point.position;
    float currentDistance = delta.length();
    
    // as our vertices haven't moved we only need to check any new vertices
    for (v = mLastNumOfVerts; v < mVertices.size(); v++) {
      delta = mVertices[v] - point.position;
      float distance = delta.length();
      if (distance < currentDistance) {
        // this one is now our closest
        point.closestVertice = v;
        currentDistance = distance;
      };
    };
    
    if (currentDistance < pCutOffDistance) {
      // we're done with this one...
      mAttractionPoints.erase(mAttractionPoints.begin() + i);
    } else {
      // copy back our new closest vertice and advance...
      mAttractionPoints[i].closestVertice = point.closestVertice;
      
      if (currentDistance < pMaxDistance) {
        // count our vertice
        numOfAPoints[point.closestVertice] += 1.0;
        vec3 norm = mAttractionPoints[i].position - mVertices[point.closestVertice];
        directions[point.closestVertice] += norm.normalize();        
      };
      
      // and advance
      i++;
    };
  };
  
  // Update our last number of vertices
  mLastNumOfVerts = numVerts;
  
  // Now check which vertices need to branch out...
  for (v = 0; v < numVerts; v++) {    
    if (numOfAPoints[v] > 0.0) {
      vec3  vert = mVertices[v];
      directions[v] /= numOfAPoints[v];
      directions[v] = directions[v].normalize() * pBranchSize;        
      vert += directions[v] + pBias;
      
      growBranch(v, vert);      
    };
  };
  
  // as long as we still have attraction points left we must still be growing our tree
  return mAttractionPoints.size() > 0; 
};

I've added a few optimisations into this code.
First of all is that I only check new vertices added since the last iteration and remember the closest vertex for the next round. As vertices once generated do not move only new vertices can be closer then the current closest vertex.
Second, instead of building a set of attraction points that are closest to a particular vertex I simply count the number of attraction points and add their normalised vectors together. Then once I'm finished I simply divide the sum of the vectors by the number of points to get the direction in which to grow my branch.

The end result is a nice stick figure tree:


It's a lot of fun experimenting with the parameters and seeing what different shapes to come up with. One of the things I added in the sample code is to generate a second point cloud after the tree has been build to create roots.

More next time....