Create a gist now

Instantly share code, notes, and snippets.

@mitjmcc /0_Code Samples.md Secret
Last active Apr 13, 2018

What would you like to do?

Collection of samples of code for effects and procedural generation

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.

oilpaint

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
}
}
}
#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();
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment