OpenCL / CUDA konsult

Söker ni efter en erfaren konsult med kunskaper i OpenCL / CUDA programmering? Kontakta mig!

Här är några exempel på vad jag jobbat med;
AI (Neurala nätverk): http://sharprbm.codeplex.com/
Shaders: https://lotsacode.wordpress.com/2013/04/17/shadertoy/
Prtikelsimulering: https://lotsacode.wordpress.com/2013/04/16/fun-with-particles-yet-another-shader-editor/

ShaderToy

ShaderToy is an site for writing and running WebGL shaders online.

Here are two screenshots from a Kaleidoscope shader I created;

image

image

 

Here’s the basic OpenCL code you’ll need to create your own kaleidoscope;

const int numPoints = 12;
const bool showFolds = false;

float rand( vec2 n ) {
	return fract(sin(dot(n.xy, vec2(12.9898, 78.233)))* 43758.5453);
}

struct Ray
{
	vec2 point;
	vec2 direction;
};

float noise(vec2 n) {
	const vec2 d = vec2(0.0, 1.0);
	vec2 b = floor(n), f = smoothstep(vec2(0.0), vec2(1.0), fract(n));
	return mix(mix(rand(b), rand(b + d.yx), f.x), mix(rand(b + d.xy), rand(b + d.yy), f.x), f.y);
}

vec2 noise2(vec2 n)
{
	return vec2(noise(vec2(n.x+0.2, n.y-0.6)), noise(vec2(n.y+3., n.x-4.)));
}

Ray GetRay(float i)
{
	vec2 position = noise2(vec2(i*6.12+iGlobalTime*0.1, i*4.43+iGlobalTime*0.1));
	return Ray(
		position,
		normalize(noise2(vec2(i*7.+iGlobalTime*0.05, i*6.))*2.0-1.0));	
}

void main(void)
{
	vec2 uv = gl_FragCoord.xy / min(iResolution.x,iResolution.y);
	vec2 curPos = uv;
	
	for(int i=0;i<numPoints;i++)
	{
		Ray ray=GetRay(float(i+1)*3.);	
			
		if(length(ray.point-curPos)<0.01 && showFolds)
		{
			gl_FragColor.rgb = vec3(1,1,1);
			return;
		}
		else if (length(curPos-(ray.point+ray.direction*0.1))<0.01 && showFolds)
		{
			gl_FragColor.rgb = vec3(1,0,0);
			return;
		}
		else
		{
			float offset=dot(curPos-ray.point, ray.direction);
			if(abs(offset)<0.001 && showFolds)
			{
				gl_FragColor.rgb = vec3(0,0,1);
				return;
			}
			if(offset<0.)
			{
				curPos -= ray.direction*offset*2.0;
			}									
		}
	}
	
	gl_FragColor.rgb = texture2D( iChannel0, vec2(curPos.x,curPos.y) ).xyz;	
	gl_FragColor.a = 1.0;		
}

Fun with Particles; Yet Another Shader Editor

During my previous experiment, with particles on a Perlin Forcefield, I came across Yet Another Shader Editor / http://yase.chnk.us/, really fun site (to me) where you can enter some OpenCL code to control a particle and a lot of particles will be rendered which follow the code. The site was created by Ryan Alexander.

On this page I intend to collect my YOSE experiments;

Particles on perlin force field:http://yase.chnk.us/#zz

image

Phlegm fountain: http://yase.chnk.us/#z7t

image

 

#define sim_res 1700
#define pixel_scale 2

const float timeVariation=0.03; const float worldSize = 10.;
const float speed=0.0055; const float spawnSize = 0.1;
const float maxAge=600.; const float outwardsForce = .2;
const float aging=0.002;

vec3 perlinForceField(vec4 pos, float variedTime){ 
  /* A time-warying Simplex force field */
  return
    vec3( 
      noise(pos.x, pos.y, pos.z + variedTime), 
      noise(pos.x+400. + variedTime, pos.y, pos.z),
      noise(pos.x+500., pos.y + variedTime, pos.z));
}

vec3 heartbeat(
  vec4 pos, float time){
	vec3 push= normalize(vec3(pos))*0.007;
  push.z=abs(push.z)*3.;
  float beat=sin(time*5.5);
  if(beat<0.)
  {
    beat*=.3;
  }
  return push*beat;
}

void stepPos(in float i, in vec4 prevPos, 
             in vec4 pos, out vec4 nextPos) {
  float variedTime = time*timeVariation;
  vec3 delta = perlinForceField(pos, variedTime);
  delta = normalize(delta)*speed;
	delta = delta + heartbeat(pos, time);
	nextPos.x = pos.x + delta.x;
	nextPos.y = pos.y + delta.y;
	nextPos.z = pos.z + delta.z;
  nextPos.w = pos.w+aging;  

  if(abs(nextPos.x)>worldSize||abs(nextPos.y)>worldSize||
     abs(nextPos.z)>worldSize||nextPos.z<-10.0||
     nextPos.w>maxAge*aging||time<=0.01)
  {
    float range=worldSize*spawnSize;
    nextPos.x = noise(pos.x, pos.y, i*3.+time)*range;
    nextPos.y = noise(pos.x, pos.y, i*3.+time*3.)*range;
    nextPos.z = -2.+noise(pos.x, pos.y, i*3.+time*5.)*range;
    nextPos.w = noise(pos.x, pos.y, i*3.+time*5.)*maxAge*aging;
  }    
}

 

Turbulent star: http://yase.chnk.us/#2r

image

#define sim_res 1500
#define pixel_scale 1

const float timeVariation=0.07; const float worldSize = 10.;
const float speed=0.0075; const float spawnSize = 0.1;
const float outwardsForce = .2;
const float radius=3.;

vec3 perlinForceField(vec4 pos, float variedTime){ 
  /* A time-warying Simplex force field */
  vec3 force=
    vec3( 
      noise(pos.x, pos.y, pos.z + variedTime), 
      noise(pos.x+400. + variedTime, pos.y, pos.z),
      noise(pos.x+500., pos.y + variedTime, pos.z));

  force = normalize(force);  

  return force;
}

vec3 sphereForce(vec3 pos){
  float len=length(pos);
	if(len<radius)
  {
		return normalize(pos)*(radius-len);
  }

	return -normalize(vec3(pos))*0.004;
}

void stepPos(in float i, in vec4 prevPos, 
             in vec4 pos, out vec4 nextPos) {
  float variedTime = time*timeVariation;
  vec3 delta = perlinForceField(pos, variedTime)*speed; 
	delta = delta + sphereForce(pos.xyz);
	nextPos.x = pos.x + delta.x;
	nextPos.y = pos.y + delta.y;
	nextPos.z = pos.z + delta.z;

	bool dead=noise(i+frame, i-frame, i+frame)<-0.99;
  if(abs(nextPos.x)>worldSize||abs(nextPos.y)>worldSize||
     abs(nextPos.z)>worldSize||nextPos.z<-10.0||
     dead||time<=0.01)
  {
    vec3 spherePos=
      normalize(
        vec3(
          noise(pos.x, pos.y, i*1.+time),
          noise(pos.x, pos.y, i*2.+time*3.),
          noise(pos.x, pos.y, i*3.+time*5.)))*radius;    

    nextPos.xyz = spherePos.xyz;
  } 

  nextPos.w = 0.1+min(1.,(length(nextPos.xyz)-radius)/(1.1));
}

 

Hairball: http://yase.chnk.us/#1n

image

Based on one of Ryan Alexanders originals, but using perlin force fields to perturb the particles.

#define sim_res 312

vec3 spring(in vec3 a, in vec3 b, in float len, in float power) {
  vec3 o = b - a;
  float m = length(o) + 0.00001;
  return (o / m) * power * (m - len);
}

void stepPos(in float i, in vec4 prevPos, in vec4 pos, out vec4 nextPos) {
	float	n = 500.;
  float t = i / count, noc = n / count;
  float g = floor(t * n), gt = fract(t * n), omgt = 1. - gt;

	int ci = int(count), ii = int(i);

  // Step
  vec3 np = pos.xyz + (pos.xyz - prevPos.xyz) * 0.9;

	// Spring
  float jf, pwr = .01, len = 1. / (count / n);
  vec3 p;
  for(int j = 1; j < 16; ++j) {
    jf = float(j);
    if(ii - j > 0 && floor((i - jf) * noc) == g) {
      p = getPos(i - jf).xyz;
    	np += spring(pos.xyz, p, len * jf, pwr);
    }
    if(ii + j < ci - 1 && floor((i + jf) * noc) == g) {
      p = getPos(i + jf).xyz;
    	np += spring(pos.xyz, p, len * jf, pwr);
    }
  }

  // Turbulence

  //nextPos.xyz = mix(rand3(g + floor(frame / 500.)) * omgt, np, 1. - gt * gt * 0.05);


  vec3 turbulence = 
    (0.5+sin(time*2.))*
    vec3(
      noise(pos.x*2.+ frame*0.005, pos.y, pos.y )*0.001,
      noise(pos.x*2., pos.y+ frame*0.005, pos.z+100.)*0.001,
      noise(pos.x*2., pos.y, pos.z+200.+ frame*0.005)*0.001);

  nextPos.xyz = mix(rand3(g) * omgt, np, 1. - gt * gt * 0.05)+turbulence;
  nextPos.w = omgt * omgt;
}

Experimental

Particles bouncing of sphere and plane: http://yase.chnk.us/#lx

Colored Perlin Particle Field

I was inspired and wanted to play with perlin particle fields, where perlin noise is used to simulate a particle moving along a force field. I found the excellent tool Processing that makes it really easy to play with stuff like this. As luck would have it, someone had already implemented Perlin Particle Fields, but in black and white.

image

I wanted color. So a few tweaks later and behold! I create particles of 7 different colors and they behave subtly different – they follow subtly different force fields. The force fields are similar enough to be fairly cohesive, but they strive for separation (see the colored “channels” at the bottom left part of the image below).

image

 

Movies

Here’s a movie of how a field like the one above is constructed, from the time particles are born until they’ve moved for a couple of thousand updates.

Colored Perlin Particle Field

 

Cooler movies (?)

But then I found Perlin noise + vectorfield + particles = fun and decided I want to have a go at that too! What he does is run his particles for a set number of iterations, save one (movie)frame, reset with a slightly different Perlin noise field, and then generate a new image – stitching them together to a movie. I’m doing the same below. I randomize the particles using the same random seed, so the particles keep getting spawned in the same locations. Without random spawn points, the movie becomes very chaotic.

All these were created using MakeAVI.

30000 articles, 1000 steps

Particles = colored points, Steps=the particles were drawn and moved.

1000 particles, 1000 steps

200 particles, 1000 steps

WebGL

Here’s a 3D WebGL version I created at Yet Another Shader Editor. You’ll need a WebGL capable graphics card to run it.

image

If you liked that one, try Flegm fountain.

Simplex Noise in C#

In a previous post I published C# code for Perlin Noise. This is C# for Simplex Noise, which is actually better and faster. Even Ken Perlin says so – and he actually invented both Simplex Noise and Perlin Noise.

 

    public class SimplexNoise
    {
        private const double Div6 = 1.0 / 6.0;
        private const double Div3 = 1.0 / 3.0;

        /* Unit vectors for gradients to points on cube,equal distances apart (ie vector from center to the middle of each side */
        private static readonly int[][] Grad3 = new[]
            {
                new[] { 1, 1, 0 }, new[] { -1, 1, 0 }, new[] { 1, -1, 0 }, new[] { -1, -1, 0 },
                new[] { 1, 0, 1 }, new[] { -1, 0, 1 }, new[] { 1, 0, -1 }, new[] { -1, 0, -1 },
                new[] { 0, 1, 1 }, new[] { 0, -1, 1 }, new[] { 0, 1, -1 }, new[] { 0, -1, -1 }
            };

        //0..255, randomized
        private static readonly int[] P = 
            {
                151, 160, 137, 91, 90, 15, 131, 13, 201, 95, 96, 53, 194, 233, 7, 225, 140, 36, 103,
                30, 69, 142, 8, 99, 37, 240, 21, 10, 23, 190, 6, 148, 247, 120, 234, 75, 0, 26, 197,
                62, 94, 252, 219, 203, 117, 35, 11, 32, 57, 177, 33, 88, 237, 149, 56, 87, 174, 20,
                125, 136, 171, 168, 68, 175, 74, 165, 71, 134, 139, 48, 27, 166, 77, 146, 158, 231,
                83, 111, 229, 122, 60, 211, 133, 230, 220, 105, 92, 41, 55, 46, 245, 40, 244, 102,
                143, 54, 65, 25, 63, 161, 1, 216, 80, 73, 209, 76, 132, 187, 208, 89, 18, 169, 200,
                196, 135, 130, 116, 188, 159, 86, 164, 100, 109, 198, 173, 186, 3, 64, 52, 217, 226,
                250, 124, 123, 5, 202, 38, 147, 118, 126, 255, 82, 85, 212, 207, 206, 59, 227, 47, 16,
                58, 17, 182, 189, 28, 42, 223, 183, 170, 213, 119, 248, 152, 2, 44, 154, 163, 70,
                221, 153, 101, 155, 167, 43, 172, 9, 129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224,
                232, 178, 185, 112, 104, 218, 246, 97, 228, 251, 34, 242, 193, 238, 210, 144, 12,
                191, 179, 162, 241, 81, 51, 145, 235, 249, 14, 239, 107, 49, 192, 214, 31, 181, 199,
                106, 157, 184, 84, 204, 176, 115, 121, 50, 45, 127, 4, 150, 254, 138, 236, 205, 93,
                222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 215, 61, 156, 180
            };

        private static readonly int[] Perm = new int[512];

        private static readonly double Sqrt3 = Math.Sqrt(3);
        private static readonly double F2 = 0.5 * (Sqrt3 - 1.0);
        private static readonly double G2 = (3.0 - Sqrt3) * Div6;

        static SimplexNoise()
        {
            for (int i = 0; i < 512; i++)
            {
                Perm[i] = P[i & 255];
            }

            Singleton = new SimplexNoise();
        }

        public static SimplexNoise Singleton { get; private set; }

        public double Noise01(double x, double y)
        {
            return (Noise(x, y) + 1) * 0.5;
        }

        public double MultiNoise(int octaves, double x, double y)
        {
            double value = 0.0;
            float mul = 1;
            for (int i = 0; i < octaves; i++)
            {
                value += Noise((x + 10) * mul, (y + 15) * mul) / mul;

                mul *= 2;
            }
            return value;
        }

        public double MultiNoise01(int octaves, double x, double y)
        {
            return (MultiNoise(octaves, x, y) + 1.0) * 0.5;
        }

        public double RidgedMulti(int octaves, double x, double y)
        {
            double value = 0.0;
            double mul = 1;
            for (int i = 0; i < octaves; i++)
            {
                double added = Noise(x * mul, y * mul) / mul;
                value += Math.Abs(added);

                mul *= 2.18387276;
            }
            return value;
        }

        public double Noise(double xin, double yin)
        {
            double n0, n1, n2; // Noise contributions from the three corners

            // Skew the input space to a square to determine which simplex cell we're in

            double s = (xin + yin) * F2; // Hairy factor for 2D
            int i = FastFloor(xin + s);
            int j = FastFloor(yin + s);

            double t = (i + j) * G2;
            double x0p = i - t; // Unskew the cell origin back to (x,y) space
            double y0p = j - t;
            double x0 = xin - x0p; // The x,y distances from the cell origin
            double y0 = yin - y0p;

            // For the 2D case, the simplex shape is an equilateral triangle.
            // Determine which simplex we are in.
            int i1, j1; // Offsets for second (middle) corner of simplex in (i,j) coords
            if (x0 > y0)
            {
                // lower triangle, XY order: (0,0)->(1,0)->(1,1)
                i1 = 1;
                j1 = 0;
            }
            else
            {
                i1 = 0;
                j1 = 1;
            } // upper triangle, YX order: (0,0)->(0,1)->(1,1)

            // A step of (1,0) in (i,j) means a step of (1-c,-c) in (x,y), and
            // a step of (0,1) in (i,j) means a step of (-c,1-c) in (x,y), where
            // c = (3-sqrt(3))/6

            double x1 = x0 - i1 + G2; // Offsets for middle corner in (x,y) unskewed coords
            double y1 = y0 - j1 + G2;
            double x2 = x0 - 1.0 + 2.0 * G2; // Offsets for last corner in (x,y) unskewed coords
            double y2 = y0 - 1.0 + 2.0 * G2;

            // Work out the hashed gradient indices of the three simplex corners
            int ii = i & 255;
            int jj = j & 255;
            int gi0 = Perm[ii + Perm[jj]] % 12;
            int gi1 = Perm[ii + i1 + Perm[jj + j1]] % 12;
            int gi2 = Perm[ii + 1 + Perm[jj + 1]] % 12;

            // Calculate the contribution from the three corners
            double t0 = 0.5 - x0 * x0 - y0 * y0;
            if (t0 < 0)
            {
                n0 = 0.0;
            }
            else
            {
                t0 *= t0;
                n0 = t0 * t0 * Dot(Grad3[gi0], x0, y0); // (x,y) of grad3 used for 2D gradient
            }
            double t1 = 0.5 - x1 * x1 - y1 * y1;
            if (t1 < 0)
            {
                n1 = 0.0;
            }
            else
            {
                t1 *= t1;
                n1 = t1 * t1 * Dot(Grad3[gi1], x1, y1);
            }
            double t2 = 0.5 - x2 * x2 - y2 * y2;
            if (t2 < 0)
            {
                n2 = 0.0;
            }
            else
            {
                t2 *= t2;
                n2 = t2 * t2 * Dot(Grad3[gi2], x2, y2);
            }
            // Add contributions from each corner to get the final noise value.
            // The result is scaled to return values in the interval [-1,1].
            return 70.0 * (n0 + n1 + n2);
        }

        public double Multi01(int octaves, double x, double y, double z)
        {
            return (Multi(octaves, x, y, z) + 1) * 0.5;
        }

        public double Multi(int octaves, double x, double y, double z)
        {
            double value = 0.0;
            double mul = 1;
            for (int i = 0; i < octaves; i++)
            {
                double added = Noise(x * mul, y * mul, z * mul) / mul;
                value += added;
                mul *= 2;
            }
            return value;
        }

        public double Noise01(double x, double y, double z)
        {
            // Noise  is in the range -1 to +1
            double val = Noise(x, y, z);
            return (val + 1) * 0.5;
        }

        public double RidgedMulti(int octaves, double x, double y, double z)
        {
            double value = 0.0;
            double mul = 1;
            for (int i = 0; i < octaves; i++)
            {
                double added = Noise(x * mul, y * mul, z * mul) / mul;
                value += Math.Abs(added);
                mul *= 2;
            }
            return value;
        }

        public double Noise(double xin, double yin, double zin)
        {
            double n0, n1, n2, n3; // Noise contributions from the four corners
            // Skew the input space to determine which simplex cell we're in

            double s = (xin + yin + zin) * Div3; // Very nice and simple skew factor for 3D
            int i = FastFloor(xin + s);
            int j = FastFloor(yin + s);
            int k = FastFloor(zin + s);

            double t = (i + j + k) * Div6;
            double ax0 = i - t; // Unskew the cell origin back to (x,y,z) space
            double ay0 = j - t;
            double az0 = k - t;
            double x0 = xin - ax0; // The x,y,z distances from the cell origin
            double y0 = yin - ay0;
            double z0 = zin - az0;
            // For the 3D case, the simplex shape is a slightly irregular tetrahedron.
            // Determine which simplex we are in.
            int i1, j1, k1; // Offsets for second corner of simplex in (i,j,k) coords
            int i2, j2, k2; // Offsets for third corner of simplex in (i,j,k) coords
            if (x0 >= y0)
            {
                if (y0 >= z0)
                {
                    i1 = 1;
                    j1 = 0;
                    k1 = 0;
                    i2 = 1;
                    j2 = 1;
                    k2 = 0;
                }
                else if (x0 >= z0)
                {
                    i1 = 1;
                    j1 = 0;
                    k1 = 0;
                    i2 = 1;
                    j2 = 0;
                    k2 = 1;
                }
                else
                {
                    i1 = 0;
                    j1 = 0;
                    k1 = 1;
                    i2 = 1;
                    j2 = 0;
                    k2 = 1;
                }
            }
            else
            {
                // x0<y0
                if (y0 < z0)
                {
                    i1 = 0;
                    j1 = 0;
                    k1 = 1;
                    i2 = 0;
                    j2 = 1;
                    k2 = 1;
                }
                else if (x0 < z0)
                {
                    i1 = 0;
                    j1 = 1;
                    k1 = 0;
                    i2 = 0;
                    j2 = 1;
                    k2 = 1;
                }
                else
                {
                    i1 = 0;
                    j1 = 1;
                    k1 = 0;
                    i2 = 1;
                    j2 = 1;
                    k2 = 0;
                }
            }
            // A step of (1,0,0) in (i,j,k) means a step of (1-c,-c,-c) in (x,y,z),
            // a step of (0,1,0) in (i,j,k) means a step of (-c,1-c,-c) in (x,y,z), and
            // a step of (0,0,1) in (i,j,k) means a step of (-c,-c,1-c) in (x,y,z), where
            // c = 1/6.
            double x1 = x0 - i1 + Div6; // Offsets for second corner in (x,y,z) coords
            double y1 = y0 - j1 + Div6;
            double z1 = z0 - k1 + Div6;
            double x2 = x0 - i2 + 2.0 * Div6; // Offsets for third corner in (x,y,z) coords
            double y2 = y0 - j2 + 2.0 * Div6;
            double z2 = z0 - k2 + 2.0 * Div6;
            double x3 = x0 - 1.0 + 3.0 * Div6; // Offsets for last corner in (x,y,z) coords
            double y3 = y0 - 1.0 + 3.0 * Div6;
            double z3 = z0 - 1.0 + 3.0 * Div6;
            // Work out the hashed gradient indices of the four simplex corners
            int ii = i & 255;
            int jj = j & 255;
            int kk = k & 255;
            int gi0 = Perm[ii + Perm[jj + Perm[kk]]] % 12;
            int gi1 = Perm[ii + i1 + Perm[jj + j1 + Perm[kk + k1]]] % 12;
            int gi2 = Perm[ii + i2 + Perm[jj + j2 + Perm[kk + k2]]] % 12;
            int gi3 = Perm[ii + 1 + Perm[jj + 1 + Perm[kk + 1]]] % 12;
            // Calculate the contribution from the four corners
            double t0 = 0.6 - x0 * x0 - y0 * y0 - z0 * z0;
            if (t0 < 0)
            {
                n0 = 0.0;
            }
            else
            {
                t0 *= t0;
                n0 = t0 * t0 * Dot(Grad3[gi0], x0, y0, z0);
            }
            double t1 = 0.6 - x1 * x1 - y1 * y1 - z1 * z1;
            if (t1 < 0)
            {
                n1 = 0.0;
            }
            else
            {
                t1 *= t1;
                n1 = t1 * t1 * Dot(Grad3[gi1], x1, y1, z1);
            }
            double t2 = 0.6 - x2 * x2 - y2 * y2 - z2 * z2;
            if (t2 < 0)
            {
                n2 = 0.0;
            }
            else
            {
                t2 *= t2;
                n2 = t2 * t2 * Dot(Grad3[gi2], x2, y2, z2);
            }
            double t3 = 0.6 - x3 * x3 - y3 * y3 - z3 * z3;
            if (t3 < 0)
            {
                n3 = 0.0;
            }
            else
            {
                t3 *= t3;
                n3 = t3 * t3 * Dot(Grad3[gi3], x3, y3, z3);
            }
            // Add contributions from each corner to get the final noise value.
            // The result is scaled to stay just inside [-1,1]
            return 32.0 * (n0 + n1 + n2 + n3);
        }

        private static int FastFloor(double x)
        {
            // This method is a *lot* faster than using (int)Math.floor(x)
            // This comment is regards to C, 
            return x > 0 ? (int)x : (int)x - 1;
        }

        private static double Dot(int[] g, double x, double y)
        {
            return g[0] * x + g[1] * y;
        }

        private static double Dot(int[] g, double x, double y, double z)
        {
            return g[0] * x + g[1] * y + g[2] * z;
        }    
    }

Catmull Clark Surface Subdivider in C#

I wanted to play with stuff like this; http://www.flickr.com/photos/junearch/sets/72157602857326020/detail/ which meant I needed a surface subdivider. Since it’s what I do for fun, I implemented the method described here: http://rosettacode.org/wiki/Catmull%E2%80%93Clark_subdivision_surface in C#. You’ll find the code below.

Undivided

Notice how the box has a hole – you can of course divide solid manifolds as well.

undivided

Once Divided

once divided

Twice divided

twice divided

Thrice divided

image

Four times divided

Dividing the box 4 times results in 10460 triangles. The divisions took 245ms– but I haven’t spent a lot of time optimizing this.

image

The code

using System;
using System.Collections.Generic;
using System.Linq;

namespace Subdivision.Core
{
    public class CatmullClarkSubdivider : ISubdivider
    {
        public Shape Subdivide(Shape shape)
        {
            // http://rosettacode.org/wiki/Catmull%E2%80%93Clark_subdivision_surface 
            Shape subdivided = new Shape();

            //for each face, a face point is created which is the average of all the points of the face.
            CreateFacePoints(shape);

            //for each edge, an edge point is created which is the average between the center of 
            //  the edge and the center of the segment made with the face points of the two adjacent faces.
            CreateEdgePoints(shape);

            //for each vertex point, its coordinates are updated from (new_coords):
            //    the old coordinates (old_coords),
            //    the average of the face points of the faces the point belongs to (avg_face_points),
            //    the average of the centers of edges the point belongs to (avg_mid_edges),
            //    how many faces a point belongs to (n), then use this formula: 
            //m1 = (n - 3) / n
            //m2 = 1 / n
            //m3 = 2 / n
            //new_coords = (m1 * old_coords)
            //           + (m2 * avg_face_points)
            //           + (m3 * avg_mid_edges)
            CreateVertexPoints(shape);

            //for a triangle face (a,b,c): 
            //   (a, edge_pointab, face_pointabc, edge_pointca)
            //   (b, edge_pointbc, face_pointabc, edge_pointab)
            //   (c, edge_pointca, face_pointabc, edge_pointbc)

            //for a quad face (a,b,c,d): 
            //   (a, edge_pointab, face_pointabcd, edge_pointda)
            //   (b, edge_pointbc, face_pointabcd, edge_pointab)
            //   (c, edge_pointcd, face_pointabcd, edge_pointbc)
            //   (d, edge_pointda, face_pointabcd, edge_pointcd)
            CreateFaces(shape, subdivided);

            return subdivided;
        }

        private void CreateFacePoints(Shape shape)
        {
            //for each face, a face point is created which is the average of all the points of the face.
            foreach (Face face in shape.Faces)
            {
                List points = face.AllPoints;
                face.FacePoint = new Point(Average(points));
            }
        }

        private void CreateEdgePoints(Shape shape)
        {
            //for each edge, an edge point is created which is the average between the center of 
            //  the edge and the center of the segment made with the face points of the two adjacent faces.
            List edges = shape.AllEdges;
            foreach (Edge edge in edges)
            {
                if (edge.IsOnBorderOfHole)
                {
                    Vector3 position =
                        Average(
                            edge.Points[0],
                            edge.Points[1]);

                    edge.EdgePoint = new Point(position);
                }
                else
                {
                    Vector3 position =
                        Average(
                            edge.Points[0],
                            edge.Points[1],
                            edge.Faces[0].FacePoint,
                            edge.Faces[1].FacePoint);

                    edge.EdgePoint = new Point(position);
                }
            }
        }

        private void CreateVertexPoints(Shape shape)
        {
            //for each vertex point, its coordinates are updated from (new_coords):
            //    the old coordinates (old_coords),
            //    the average of the centers of edges the point belongs to (avg_mid_edges),
            //    the average of the face points of the faces the point belongs to (avg_face_points),
            //    how many faces a point belongs to (n), then use this formula: 
            //m1 = (n - 3) / n
            //m2 = 1 / n
            //m3 = 2 / n
            //new_coords = (m1 * old_coords)
            //           + (m2 * avg_face_points)
            //           + (m3 * avg_mid_edges) 
            List allPoints = shape.AllPoints;
            List allEdges = shape.AllEdges;

            foreach (Point oldPoint in allPoints)
            {
                if (oldPoint.IsOnBorderOfHole)
                {
                    oldPoint.Successor = CreateVertexPointForBorderPoint(oldPoint);
                }
                else
                {
                    oldPoint.Successor = CreateVertexPoint(allEdges, oldPoint);
                }
            }
        }

        private Point CreateVertexPoint(List allEdges, Point oldPoint)
        {
            //    the average of the face points of the faces the point belongs to (avg_face_points),
            //    how many faces a point belongs to (n), then use this formula: 

            //    the average of the centers of edges the point belongs to (avg_mid_edges),
            Vector3 avgMidEdges = Vector3.Average(oldPoint.Edges.Select(e => e.Middle));

            //    the average of the face points of the faces the point belongs to (avg_face_points),
            List pointFaces = oldPoint.Edges.SelectMany(e => e.Faces).Distinct().ToList();
            Vector3 avgFacePoints = Average(pointFaces.Select(pf => pf.FacePoint));

            int faceCount = pointFaces.Count;

            double m1 = (faceCount - 3f) / faceCount;
            double m2 = 1f / faceCount;
            double m3 = 2f / faceCount;

            Point newPoint = new Point(m1 * oldPoint.Position + m2 * avgFacePoints + m3 * avgMidEdges);
            return newPoint;
        }

        private Point CreateVertexPointForBorderPoint(Point oldPoint)
        {
            //for the vertex points that are on the border of a hole, the new coordinates are calculated as follows:
            // in all the edges the point belongs to, only take in account the middles of the edges that are on the border of the hole
            // calculate the average between these points (on the hole boundary) and the old coordinates (also on the hole boundary). 
            List positions = oldPoint.Edges.Where(e => e.IsOnBorderOfHole).Select(e => e.Middle).ToList();
            positions.Add(oldPoint.Position);

            return new Point(Vector3.Average(positions));
        }

        private void CreateFaces(Shape shape, Shape subdivided)
        {
            List faces = shape.Faces;
            List existingEdges = new List();
            foreach (Face face in faces)
            {
                if (face.AllPoints.Count() == 3)
                {
                    CreateTriangleFace(existingEdges, subdivided, face);
                }
                else if (face.AllPoints.Count() == 4)
                {
                    CreateQuadFace(existingEdges, subdivided, face);
                }
                else
                {
                    throw new InvalidOperationException(string.Format("Unhandled facetype (point count={0})!", face.AllPoints.Count()));
                }
            }
            SubdivisionUtilities.VerifyThatThereAreNoEdgeDuplicates(existingEdges);
        }

        private void CreateTriangleFace(List existingEdges, Shape subdivided, Face face)
        {
            List points = face.AllPoints;
            Point a = points[0];
            Point b = points[1];
            Point c = points[2];

            //for a triangle face (a,b,c): 
            //   (a, edge_pointab, face_pointabc, edge_pointca)
            //   (b, edge_pointbc, face_pointabc, edge_pointab)
            //   (c, edge_pointca, face_pointabc, edge_pointbc)
            Point facePoint = face.FacePoint;

            subdivided.Faces.Add(SubdivisionUtilities.CreateFaceF(existingEdges, a, face.Edges[0].EdgePoint, facePoint, face.Edges[2].EdgePoint));
            subdivided.Faces.Add(SubdivisionUtilities.CreateFaceF(existingEdges, b, face.Edges[1].EdgePoint, facePoint, face.Edges[0].EdgePoint));
            subdivided.Faces.Add(SubdivisionUtilities.CreateFaceF(existingEdges, c, face.Edges[2].EdgePoint, facePoint, face.Edges[1].EdgePoint));

            SubdivisionUtilities.VerifyThatThereAreNoEdgeDuplicates(existingEdges);
        }

        private void CreateQuadFace(List existingEdges, Shape subdivided, Face face)
        {
            //                  0 1 2 -> 3 
            //for a quad face (a,b,c,d): 
            //   (a, edge_pointab, face_pointabcd, edge_pointda)
            //   (b, edge_pointbc, face_pointabcd, edge_pointab)
            //   (c, edge_pointcd, face_pointabcd, edge_pointbc)
            //   (d, edge_pointda, face_pointabcd, edge_pointcd)
            List points = face.AllPoints;
            Point a = points[0].Successor;
            Point b = points[1].Successor;
            Point c = points[2].Successor;
            Point d = points[3].Successor;

            Point facePoint = face.FacePoint;

            subdivided.Faces.Add(SubdivisionUtilities.CreateFaceF(existingEdges, a, face.Edges[0].EdgePoint, facePoint, face.Edges[3].EdgePoint));
            subdivided.Faces.Add(SubdivisionUtilities.CreateFaceF(existingEdges, b, face.Edges[1].EdgePoint, facePoint, face.Edges[0].EdgePoint));
            subdivided.Faces.Add(SubdivisionUtilities.CreateFaceF(existingEdges, c, face.Edges[2].EdgePoint, facePoint, face.Edges[1].EdgePoint));
            subdivided.Faces.Add(SubdivisionUtilities.CreateFaceF(existingEdges, d, face.Edges[3].EdgePoint, facePoint, face.Edges[2].EdgePoint));

            SubdivisionUtilities.VerifyThatThereAreNoEdgeDuplicates(existingEdges);
        }

        private Vector3 Average(IEnumerable points)
        {
            return Vector3.Average(points.Select(p => p.Position));
        }

        private Vector3 Average(params Point[] points)
        {
            return Vector3.Average(points.Select(p => p.Position));
        }
    }
}

You can download the complete source code here (be warned, it’s 17MB).