0_Code Samples.md

Collection of samples of code for effects and procedural generation

1_Kuwahara.md

Implementation of a Kuwahara Filter Post Processing Effect

Anisotropic Kuwahara Filtering on the GPU

Jan Eric Kyprianidis, Henry Kang, Jürgen Döllner

The effect emulates oil paint with a post processing effect using Shaderlab in Unity3D.

1_Kuwahara.shader

Shader "Custom/Kuwahara"
{
    Properties
    {
        _MainTex("Texture", 2D) = "white" {}
        _Radius ("Size of the Strokes", Range(0, 10)) = 3
    }
    SubShader
    {
        Cull Off ZWrite Off ZTest Always

        Pass
        {
            CGPROGRAM
            #pragma vertex vert
            #pragma fragment frag

            #include "UnityCG.cginc"

            uniform sampler2D _MainTex;
            uniform int _Radius;

            struct v2f {
                float4 position : POSITION;
                float4 uv : TEXCOORD0;
                float4 screenPos : TEXCOORD1;
            };

            v2f_img vert(appdata_base v) {
                v2f_img o;
                o.pos = UnityObjectToClipPos(v.vertex);
                o.uv = v.texcoord;
                return o;
            }

            struct region {
                int x1, y1, x2, y2;
            };

            float4 _MainTex_TexelSize;

            float4 frag (v2f i) : SV_Target
            {
                float w = i.screenPos.w;
                float2 uv = i.uv;
                float n = float((_Radius + 1) * (_Radius + 1));
                float4 col = tex2D(_MainTex, uv);

                float3 m[4];
                float3 s[4];

                for (int k = 0; k < 4; ++k) {
                    m[k] = float3(0, 0, 0);
                    s[k] = float3(0, 0, 0);
                }

                region R[4] = {
                    {-_Radius, -_Radius,       0,       0},
                    {       0, -_Radius, _Radius,       0},
                    {       0,        0, _Radius, _Radius},
                    {-_Radius,        0,       0, _Radius}
                };

                for (int k = 0; k < 4; ++k) {
                    for (int j = R[k].y1; j <= R[k].y2; ++j) {
                        for (int i = R[k].x1; i <= R[k].x2; ++i) {
                            float3 c = tex2D(_MainTex, uv + (float2(i * _MainTex_TexelSize.x, j * _MainTex_TexelSize.y))).rgb;
                            m[k] += c;
                            s[k] += c * c;
                        }
                    }
                }

                float min = 1e+2;
                float s2;
                for (int k = 0; k < 4; ++k) {
                    m[k] /= n;
                    s[k] = abs(s[k] / n - m[k] * m[k]);

                    s2 = s[k].r + s[k].g + s[k].b;
                    if (s2 < min) {
                        min = s2;
                        col.rgb = m[k].rgb;
                    }
                }
                return col;
            }
            ENDCG
        }
    }
}

2_ConstrainedParticles.md

Implementation of Constrained Particles

Physically Based Modeling: Constrained Dynamics

Andrew Witkin, Pixar Animation Studios

2_simulator.cpp

#include "simulator.h"
#include <fstream>
#include <iostream>

using namespace std;

double _radius = 0.2;
double _distance = 0.1;

/**
 * Implementation of Physically Based Modeling Constrained Dynamics
 * This class is instantiated and run in a main.cpp using GLUT
 */
Simulator::Simulator() {
    // initialize the particles
    mParticles.resize(2);

    // Init particle positions (default is 0, 0, 0)
    mParticles[0].mPosition[0] = _radius;
    mParticles[1].mPosition[0] = _radius;
    mParticles[1].mPosition[1] = _radius - _distance;
    
    mTimeStep = 0.0003;
}

int Simulator::getNumParticles() {
    return mParticles.size();
}

Particle* Simulator::getParticle(int index) {
    return &mParticles[index];
}

double Simulator::getTimeStep() {
    return mTimeStep;
}

void Simulator::reset() {
    mParticles[0].mPosition[1] = 0.0;
    mParticles[0].mPosition[0] = 0.2;
    mParticles[1].mPosition[0] = 0.2;
    mParticles[1].mPosition[1] = -0.1;
    
    for (int i = 0; i < getNumParticles(); i++) {
        mParticles[i].mVelocity.setZero();
        mParticles[i].mAccumulatedForce.setZero();
    }
    
}

void Simulator::simulate() {
    // Applying gravity, the only external force in this system
    for (int i = 0; i < (int) mParticles.size(); i++) {
        mParticles[i].mAccumulatedForce[1] -= 9.8 * mParticles[i].mMass;
    }

    // Caching the position and velocity of particle 1 and 2
    Eigen::Vector3d x1 = mParticles[0].mPosition;
    Eigen::Vector3d x2 = mParticles[1].mPosition;
    Eigen::Vector3d v1 = mParticles[0].mVelocity;
    Eigen::Vector3d v2 = mParticles[1].mVelocity;

    // q dot, the velocity vector
    Eigen::Matrix<double, 6, 1> qd;
    qd << v1, v2; 
    
    // Q, the global force vector
    Eigen::Matrix<double, 6, 1> Q;
    Q << mParticles[0].mAccumulatedForce,
         mParticles[1].mAccumulatedForce;

    double m0 = mParticles[0].mMass,
           m1 = mParticles[1].mMass;
    // M, the mass matrix
    M.diagonal() <<  m0, m0, m0, m1, m1, m1;
    // W, the inverse mass matrix
    W = M.inverse();

    Eigen::Vector3d x1_x2 = x1 - x2;
    // J, the Jacobian matrix of the constraints C
    J <<    x1.transpose(),   0,     0,     0,
         x1_x2.transpose(), -x1_x2.transpose();

    Eigen::Vector3d v1_v2 = v1 - v2;
    // J dot, the time derivative of the Jacobian
    Jd << v1.transpose(),            0, 0, 0,
       v1_v2.transpose(), -v1_v2.transpose();
    
    // C, the constraint vector function
    Eigen:: Matrix<double, 2, 1> C;
    C << 0.5 * (x1.dot(x1) - _radius * _radius),
         0.5 * ((x1 - x2).dot(x1 - x2) - _distance * _distance);

    // C dot, the derivative constraint vector function
    Eigen:: Matrix<double, 2, 1> Cd;
    Cd << x1.dot(v1),
          x1.dot(v1) - x1.dot(v2) - x2.dot(v1) + x2.dot(v2);

    // Computing λ, the Lagrange multiplier vector
    // given JWJ^Tλ = −J˙q˙ − JWQ − ksC − kdC˙

    // Feedback terms to account for numerical drift
    double ks = 0.2,
           kd = 0.2;
    // Collecting terms from the above equation
    Eigen::Matrix<double, 2, 1> B =  (-Jd) * qd - J * W * Q - ks * C - kd * Cd;
    Eigen::Matrix<double, 2, 2> A = J * W * J.transpose();
    // Solve for λ
    Eigen::Matrix<double, 2, 1> lambda = A.inverse() * B;
    // Compute the constraint force
    Eigen::Matrix<double, 6, 1> Qhat = J.transpose() * lambda;
    // Apply the contstraint force
    mParticles[0].mAccumulatedForce += Eigen::Vector3d(Qhat(0, 0), Qhat(1, 0), Qhat(2, 0));
    mParticles[1].mAccumulatedForce += Eigen::Vector3d(Qhat(3, 0), Qhat(4, 0), Qhat(5, 0));
    // Use Semi Implicit Euler integration
    for (int i = 0; i < (int) mParticles.size(); i++) {
        mParticles[i].mVelocity += mParticles[i].mAccumulatedForce / mParticles[i].mMass * mTimeStep;
        mParticles[i].mPosition += mParticles[i].mVelocity * mTimeStep;
    }

    for (int i = 0; i < (int) mParticles.size(); i++) {
        mParticles[i].mAccumulatedForce.setZero();
    }
}

3_CatmullClark.md

Implementation of Catmull-Clark Subdivision Surface

Catmull-Clark Subdivision Surface

Edwin Catmull, Jim Clark