Chapter 02

物理ベースレンダリングの実装

mebiusbox
mebiusbox
2021.02.23に更新
このチャプターの目次
  1. 📌 Three.js について
  2. 📌 ソースコード
  3. 📌 頂点シェーダ
  4. 📌 フラグメントシェーダ
  5. 📌 サンプル
  6. 📌 最後に
  7. 📌 参考文献

WebGLのJSライブラリ「Three.js」を使って，実際に物理ベースレンダリングを行います．
Three.js には物理ベースのマテリアルが標準で用意されていますが，ここではシェーダを独自に作成してレンダリングします．
シェーダ言語は GLSL です．実装編の読者は GLSL もしくは何らかのシェーダ言語でコーディングしたことがある人を想定しています．また，Three.js について説明を行いませんので，興味のある方は調べてみてください．

📌 Three.js について

Three.js は JavaScript 3D ライブラリです．WebGL に対応しており，機能も充実しているため手軽に3Dシーンをレンダリングすることが出来ます．ライブラリは以下から入手することが出来ます．

https://threejs.org/

「download」をクリックすればサンプルコードを含めたソースコードをダウンロード出来ます．

📌 ソースコード

実装するサンプルのソースコードは以下から入手することができます．

https://github.com/mebiusbox/pbr

今回のファイル構成は次のようになっています．

pbr/
  js/
    three.min.js
    Detector.js
    loadFiles.js
    libs/dat.gui.min.js
    libs/stats.min.js
    geometries/TeapotBufferGeometry.js
    controls/OrbitControls.js
  shaders/
    pbr_vert.glsl
    pbr_frag.glsl
  index.html

js フォルダの loadFiles.js は独自のスクリプトで，それ以外は Three.js ライブラリの一部です．

📌 頂点シェーダ

まず大事なことですが，ライティング計算はカメラ座標系で行います．これは Three.js もライティング計算はカメラ座標系になっているからです．Three.js で独自のシェーダを使う場合は THREE.ShaderMaterial を使います．THREE.ShaderMaterial を使うと，カメラやオブジェクトから変換行列を計算して自動で渡してくれます．また，ジオメトリの頂点属性も位置や法線，UV値を自動で追加してくれます．これは大変便利です．もし，接法線や従法線を頂点に追加した場合は自分で attribute を宣言する必要があります．また，自分で定義したい場合は THREE.RawShaderMaterial を使います．

以下は THREE.ShaderMaterial で作ると頂点シェーダの先頭に自動で付与されるコードです．

precision highp float;
precision highp int;
#define SHADER_NAME ShaderMaterial
#define VERTEX_TEXTURES
#define GAMMA_FACTOR 2
#define MAX_BONES 0
#define BONE_TEXTURE
#define NUM_CLIPPING_PLANES 0
uniform mat4 modelMatrix;
uniform mat4 modelViewMatrix;
uniform mat4 projectionMatrix;
uniform mat4 viewMatrix;
uniform mat3 normalMatrix;
uniform vec3 cameraPosition;
attribute vec3 position;
attribute vec3 normal;
attribute vec2 uv;
#ifdef USE_COLOR
	attribute vec3 color;
#endif
#ifdef USE_MORPHTARGETS
	attribute vec3 morphTarget0;
	attribute vec3 morphTarget1;
	attribute vec3 morphTarget2;
	attribute vec3 morphTarget3;
#ifdef USE_MORPHNORMALS
	attribute vec3 morphNormal0;
	attribute vec3 morphNormal1;
	attribute vec3 morphNormal2;
	attribute vec3 morphNormal3;
#else
	attribute vec3 morphTarget4;
	attribute vec3 morphTarget5;
	attribute vec3 morphTarget6;
	attribute vec3 morphTarget7;
#endif
#endif
#ifdef USE_SKINNING
	attribute vec4 skinIndex;
	attribute vec4 skinWeight;
#endif

そして pbr_vert.glsl の内容は次のようになっています．

varying vec3 vViewPosition;
varying vec3 vNormal;
void main() {
  vec4 mvPosition = modelViewMatrix * vec4(position, 1.0);
  gl_Position = projectionMatrix * mvPosition;
  vViewPosition = -mvPosition.xyz;
  vNormal = normalMatrix * normal;
}

varying はフラグメントシェーダに渡す変数を定義します．vViewPosition は頂点位置から視点位置までのベクトルです．mvPosition は頂点座標をカメラ座標系に変換したもので，その値をマイナスにしたものを vViewPosition に入れています．カメラ座標系で原点はカメラ位置になりますので，頂点位置の符号を反転すると，頂点位置から原点へ向かうベクトル，つまり頂点位置からカメラに向かう視線ベクトルになります．
vNormal は normalMatrix で法線を変換したものです．normalMatrix は modelViewMatrix の逆行列の転置行列が入っており，normalMatrix = transpose(inverse(modelViewMatrix)) と同じです．なぜ法線を modelViewMatrix で変換してはいけないのかはここでは説明しませんが，modelViewMatrix で変換すると正しい法線にならないからです．

📌 フラグメントシェーダ

次にフラグメントシェーダです．ここで物理ベースレンダリングの処理を行っています．
まずは，Three.js がフラグメントシェーダに自動で付与するコードは次のようになっています．

precision highp float;
precision highp int;
#define SHADER_NAME ShaderMaterial
#define GAMMA_FACTOR 2
#define NUM_CLIPPING_PLANES 0
#define UNION_CLIPPING_PLANES 0
uniform mat4 viewMatrix;
uniform vec3 cameraPosition;
#define TONE_MAPPING
#define saturate(a) clamp( a, 0.0, 1.0 )
uniform float toneMappingExposure;
uniform float toneMappingWhitePoint;
vec3 LinearToneMapping( vec3 color ) {
	return toneMappingExposure * color;
}
vec3 ReinhardToneMapping( vec3 color ) {
	color *= toneMappingExposure;
	return saturate( color / ( vec3( 1.0 ) + color ) );
}
#define Uncharted2Helper( x ) max( ( ( x * ( 0.15 * x + 0.10 * 0.50 ) + 0.20 * 0.02 ) / ( x * ( 0.15 * x + 0.50 ) + 0.20 * 0.30 ) ) - 0.02 / 0.30, vec3( 0.0 ) )
vec3 Uncharted2ToneMapping( vec3 color ) {
	color *= toneMappingExposure;
	return saturate( Uncharted2Helper( color ) / Uncharted2Helper( vec3( toneMappingWhitePoint ) ) );
}
vec3 OptimizedCineonToneMapping( vec3 color ) {
	color *= toneMappingExposure;
	color = max( vec3( 0.0 ), color - 0.004 );
	return pow( ( color * ( 6.2 * color + 0.5 ) ) / ( color * ( 6.2 * color + 1.7 ) + 0.06 ), vec3( 2.2 ) );
}

vec3 toneMapping( vec3 color ) { return LinearToneMapping( color ); }

vec4 LinearToLinear( in vec4 value ) {
	return value;
}
vec4 GammaToLinear( in vec4 value, in float gammaFactor ) {
	return vec4( pow( value.xyz, vec3( gammaFactor ) ), value.w );
}
vec4 LinearToGamma( in vec4 value, in float gammaFactor ) {
	return vec4( pow( value.xyz, vec3( 1.0 / gammaFactor ) ), value.w );
}
vec4 sRGBToLinear( in vec4 value ) {
	return vec4( mix( pow( value.rgb * 0.9478672986 + vec3( 0.0521327014 ), vec3( 2.4 ) ), value.rgb * 0.0773993808, vec3( lessThanEqual( value.rgb, vec3( 0.04045 ) ) ) ), value.w );
}
vec4 LinearTosRGB( in vec4 value ) {
	return vec4( mix( pow( value.rgb, vec3( 0.41666 ) ) * 1.055 - vec3( 0.055 ), value.rgb * 12.92, vec3( lessThanEqual( value.rgb, vec3( 0.0031308 ) ) ) ), value.w );
}
vec4 RGBEToLinear( in vec4 value ) {
	return vec4( value.rgb * exp2( value.a * 255.0 - 128.0 ), 1.0 );
}
vec4 LinearToRGBE( in vec4 value ) {
	float maxComponent = max( max( value.r, value.g ), value.b );
	float fExp = clamp( ceil( log2( maxComponent ) ), -128.0, 127.0 );
	return vec4( value.rgb / exp2( fExp ), ( fExp + 128.0 ) / 255.0 );
}
vec4 RGBMToLinear( in vec4 value, in float maxRange ) {
	return vec4( value.xyz * value.w * maxRange, 1.0 );
}
vec4 LinearToRGBM( in vec4 value, in float maxRange ) {
	float maxRGB = max( value.x, max( value.g, value.b ) );
	float M      = clamp( maxRGB / maxRange, 0.0, 1.0 );
	M            = ceil( M * 255.0 ) / 255.0;
	return vec4( value.rgb / ( M * maxRange ), M );
}
vec4 RGBDToLinear( in vec4 value, in float maxRange ) {
	return vec4( value.rgb * ( ( maxRange / 255.0 ) / value.a ), 1.0 );
}
vec4 LinearToRGBD( in vec4 value, in float maxRange ) {
	float maxRGB = max( value.x, max( value.g, value.b ) );
	float D      = max( maxRange / maxRGB, 1.0 );
	D            = min( floor( D ) / 255.0, 1.0 );
	return vec4( value.rgb * ( D * ( 255.0 / maxRange ) ), D );
}
const mat3 cLogLuvM = mat3( 0.2209, 0.3390, 0.4184, 0.1138, 0.6780, 0.7319, 0.0102, 0.1130, 0.2969 );
vec4 LinearToLogLuv( in vec4 value )  {
	vec3 Xp_Y_XYZp = value.rgb * cLogLuvM;
	Xp_Y_XYZp = max(Xp_Y_XYZp, vec3(1e-6, 1e-6, 1e-6));
	vec4 vResult;
	vResult.xy = Xp_Y_XYZp.xy / Xp_Y_XYZp.z;
	float Le = 2.0 * log2(Xp_Y_XYZp.y) + 127.0;
	vResult.w = fract(Le);
	vResult.z = (Le - (floor(vResult.w*255.0))/255.0)/255.0;
	return vResult;
}
const mat3 cLogLuvInverseM = mat3( 6.0014, -2.7008, -1.7996, -1.3320, 3.1029, -5.7721, 0.3008, -1.0882, 5.6268 );
vec4 LogLuvToLinear( in vec4 value ) {
	float Le = value.z * 255.0 + value.w;
	vec3 Xp_Y_XYZp;
	Xp_Y_XYZp.y = exp2((Le - 127.0) / 2.0);
	Xp_Y_XYZp.z = Xp_Y_XYZp.y / value.y;
	Xp_Y_XYZp.x = value.x * Xp_Y_XYZp.z;
	vec3 vRGB = Xp_Y_XYZp.rgb * cLogLuvInverseM;
	return vec4( max(vRGB, 0.0), 1.0 );
}

vec4 mapTexelToLinear( vec4 value ) { return LinearToLinear( value ); }
vec4 envMapTexelToLinear( vec4 value ) { return LinearToLinear( value ); }
vec4 emissiveMapTexelToLinear( vec4 value ) { return LinearToLinear( value ); }
vec4 linearToOutputTexel( vec4 value ) { return LinearToLinear( value ); }

ビュー行列，カメラ位置，トーンマッピング，ガンマ変換，テクスチャの RGBE, RGBM, RGBD 形式の変換などが定義されています．

幾何情報

ライティングの計算に必要なものは，「物体表面の位置と法線，視線への向き」「放射照度」「物体表面の材質」です．その内「物体表面の位置と法線，視線への向き」は幾何情報のことで，次のように定義します．

struct GeometricContext {
  vec3 position;
  vec3 normal;
  vec3 viewDir;
};

position は位置，normal は法線，viewDir は視線への向きです．
頂点シェーダから渡される vViewPosition と vNormal で幾何情報を設定します．

GeometricContext geometry;
geometry.position = -vViewPosition;
geometry.normal = normalize(vNormal);
geometry.viewDir = normalize(vViewPosition);

vViewPosition はカメラ座標での，表面位置からカメラへの視線ベクトルでした．それをまたマイナスにすると表面位置になります．そして，視線ベクトルと法線は正規化しておきます．

放射照度

レンダリング方程式では放射照度（放射輝度とコサイン項の積）が必要です．今回扱う光源は平行光源，点光源，スポットライトです．この光源からある点（＝１ピクセル）に入ってくる放射輝度を求めます．平行光源からの入射光はどの位置からでも入射してくる方向は同じですが，点光源とスポットライトは光源の位置から入射してきます．これらの光源は物体表面に光源の強さ分の放射輝度が放射されるように設定しています．ですが，例えば拡散反射のランバートモデルで考えてみると，\pi で割られるので，光源の強さよりも弱い放射輝度が入射することになります．結論として，これらの光源に対して，\pi を掛けることで調整します．この平行光源や点光源，スポットライトなどのことを punctual light と言います．

入射光（ここでは光源から入ってくる光）は IncidentLight と呼びました．入射光は入射してくる向き，光の波長（＝色）の情報を持っています．これを定義すると

struct IncidentLight {
  vec3 direction;
  vec3 color;
};

となります．
平行光源，点光源，スポットライトから入射光を求められれば，入射光から放射照度の計算とを分けることができます．点光源とスポットライトは減衰しますので，光源の位置と物体表面の位置から減衰率を求める必要があります．また，光源から光が届かないとわかれば，その光源のライティング計算をする必要がなくなりますので，IncidentLight に visible 項目を追加し，光が届くときに真となるようにします．

まず，平行光源，点光源，スポットライトの定義をします．

struct DirectionalLight {
  vec3 direction;
  vec3 color;
};

struct PointLight {
  vec3 position;
  vec3 color;
  float distance;
  float decay;
};

struct SpotLight {
  vec3 position;
  vec3 direction;
  vec3 color;
  float distance;
  float decay;
  float coneCos;
  float penumbraCos;
};

position は光源の位置，color は光の色，distance は光源からの光が届く距離，decay は減衰率，coneCos は光源から出射する光線の幅，penumbraCos は光源から出射する光線の減衰幅です．coneCos と penumbraCos はシェーダに渡すときに Math.cos で計算したものになっています．

各光源から入射光を求めるコードは次のようになっています．

bool testLightInRange(const in float lightDistance, const in float cutoffDistance) {
  return any(bvec2(cutoffDistance == 0.0, lightDistance < cutoffDistance));
}

float punctualLightIntensityToIrradianceFactor(const in float lightDistance, const in float cutoffDistance, const in float decayExponent) {
  if (decayExponent > 0.0) {
    return pow(saturate(-lightDistance / cutoffDistance + 1.0), decayExponent);
  }
  
  return 1.0;
}

void getDirectionalDirectLightIrradiance(const in DirectionalLight directionalLight, const in GeometricContext geometry, out IncidentLight directLight) {
  directLight.color = directionalLight.color;
  directLight.direction = directionalLight.direction;
  directLight.visible = true;
}

void getPointDirectLightIrradiance(const in PointLight pointLight, const in GeometricContext geometry, out IncidentLight directLight) {
  vec3 L = pointLight.position - geometry.position;
  directLight.direction = normalize(L);
  
  float lightDistance = length(L);
  if (testLightInRange(lightDistance, pointLight.distance)) {
    directLight.color = pointLight.color;
    directLight.color *= punctualLightIntensityToIrradianceFactor(lightDistance, pointLight.distance, pointLight.decay);
    directLight.visible = true;
  } else {
    directLight.color = vec3(0.0);
    directLight.visible = false;
  }
}

void getSpotDirectLightIrradiance(const in SpotLight spotLight, const in GeometricContext geometry, out IncidentLight directLight) {
  vec3 L = spotLight.position - geometry.position;
  directLight.direction = normalize(L);
  
  float lightDistance = length(L);
  float angleCos = dot(directLight.direction, spotLight.direction);
  
  if (all(bvec2(angleCos > spotLight.coneCos, testLightInRange(lightDistance, spotLight.distance)))) {
    float spotEffect = smoothstep(spotLight.coneCos, spotLight.penumbraCos, angleCos);
    directLight.color = spotLight.color;
    directLight.color *= spotEffect * punctualLightIntensityToIrradianceFactor(lightDistance, spotLight.distance, spotLight.decay);
    directLight.visible = true;
  } else {
    directLight.color = vec3(0.0);
    directLight.visible = false;
  }
}

光源からの光が届くかどうかは testLightInRange で，光源からの減衰率計算は punctualLightIntensityToIrradianceFactor で行っています．
各光源から求めた入射光に，\pi を掛けて，さらにコサイン項を掛けると放射照度になります．

物体表面の材質

物体表面の材質は「入射光のうち拡散反射になる割合＝ディフューズリフレクタンス」「入射光のうち鏡面反射になる割合＝スペキュラーリフレクタンス」「物体表面の粗さ」になります．これを定義すると次のようになります．

struct Material {
  vec3 diffuseColor;
  vec3 specularColor;
  float specularRoughness;
};

diffuseColor はディフューズリフレクタンス，specularColor はスペキュラーリフレクタンス，specularRoughness は表面の粗さです．

これらは，metallic，roughness，albedo というパラメータから計算します．metallic は金属か非金属かを直感的に設定するパラメータで，鏡面反射率です．roughness は物体表面の粗さです．albedo は波長ごとの反射能で，つまりは色（RGB）ごとの反射能です．この当たりのコードは次のようになっています．

Material material;
material.diffuseColor = mix(albedo, vec3(0.0), metallic);
material.specularColor = mix(vec3(0.04), albedo, metallic);
material.specularRoughness = roughness;

スペキュラーリフレクタンスの 0.04 は，非金属でも 4%(0.04) は鏡面反射させるということです．4% という数字は一般的な不導体を網羅する値だそうです．

BRDF

次に拡散反射BRDFと鏡面反射BRDFの実装です．ここでは前回説明した BRDF を使います．つまり，拡散反射BRDFには正規化ランバートモデルを，鏡面反射BRDFにはクック・トランスのマイクロファセットモデルを使います．マイクロファセットモデルの分布関数は GGX, 幾何減衰項は Smith モデルの Schlick 近似，フレネルには Schlick の近似式を使います．

各数式をまとめておくと

拡散反射BRDF（正規化ランバートモデル）

f_{d} = \frac{\rho_{d}}{\pi}

鏡面反射BRDF（クック・トランスのマイクロファセットモデル）

f_{r} = \frac{D(h) G(l,v) F(v,h)}{4(n \cdot l)(n \cdot v)}

マイクロファセット分布関数

D_{GGX}(h) = \frac{\alpha^2}{\pi((n \cdot h)^2(\alpha^2-1)+1)^2}

幾何減衰項

\begin{aligned} k &= \frac{\alpha}{2} \\[2ex] G_{Schlick}(v) &= \frac{n \cdot v}{(n \cdot v)(1-k) + k} \\[2ex] G_{Smith}(l,v) &= G_{Schlick}(l) G_{Schlick}(v) \end{aligned}

フレネル

F_{Schlick}(v,h) = F_{0} + (1-F_{0})(1- v \cdot h)^5

これらをコード化すると次のようになります．

// Normalized Lambert
vec3 DiffuseBRDF(vec3 diffuseColor) {
  return diffuseColor / PI;
}

vec3 F_Schlick(vec3 specularColor, vec3 H, vec3 V) {
  return (specularColor + (1.0 - specularColor) * pow(1.0 - saturate(dot(V,H)), 5.0));
}

float D_GGX(float a, float dotNH) {
  float a2 = a*a;
  float dotNH2 = dotNH*dotNH;
  float d = dotNH2 * (a2 - 1.0) + 1.0;
  return a2 / (PI * d * d);
}

float G_Smith_Schlick_GGX(float a, float dotNV, float dotNL) {
  float k = a*a*0.5 + EPSILON;
  float gl = dotNL / (dotNL * (1.0 - k) + k);
  float gv = dotNV / (dotNV * (1.0 - k) + k);
  return gl*gv;
}

// Cook-Torrance
vec3 SpecularBRDF(const in IncidentLight directLight, const in GeometricContext geometry, vec3 specularColor, float roughnessFactor) {
  
  vec3 N = geometry.normal;
  vec3 V = geometry.viewDir;
  vec3 L = directLight.direction;
  
  float dotNL = saturate(dot(N,L));
  float dotNV = saturate(dot(N,V));
  vec3 H = normalize(L+V);
  float dotNH = saturate(dot(N,H));
  float dotVH = saturate(dot(V,H));
  float dotLV = saturate(dot(L,V));
  float a = roughnessFactor * roughnessFactor;

  float D = D_GGX(a, dotNH);
  float G = G_Smith_Schlick_GGX(a, dotNV, dotNL);
  vec3 F = F_Schlick(specularColor, V, H);
  return (F*(G*D))/(4.0*dotNL*dotNV+EPSILON);
}

所々に EPSILON があります．これはゼロ除算を防ぐためのおまじないと思ってください．

レンダリング方程式

点 x から \vec{\omega} 方向に出射される（outgoing）放射輝度 L_{o} は，放射される（emitted）放射輝度 L_{e} と，反射される（reflected）放射輝度 L_{r} の和でした．

L_{o}(x,\vec{\omega}) = L_{e}(x,\vec{\omega}) + L_{r}(x,\vec{\omega})

そして，レンダリング方程式は以下のとおりです．

L_{o}(x,\vec{\omega}) = L_{e}(x,\vec{\omega}) + \int_{\Omega} f_{r}(x,\vec{\omega}',\vec{\omega}) L_{i}(x,\vec{\omega}') (\vec{\omega}' \cdot \vec{n}) d\vec{\omega}

レンダリング方程式の反射成分は，BRDF，放射輝度，コサイン項です．コサイン項は物体表面の法線とライトベクトルの内積なので，幾何情報と入射光から求まります．これで必要な情報は全て揃いました．

tuto-pbr-rendering-equation-reflection

反射成分は光源の数だけ求めることになります．反射される放射輝度を ReflectedLight として定義します．

struct ReflectedLight {
  vec3 directDiffuse;
  vec3 directSpecular;
};

directDiffuse は拡散反射成分，directSpecular は鏡面反射成分です．また，将来的な話になりますがグローバルイルミネーションも考慮して，間接光による拡散反射成分と鏡面反射成分も追加します．

struct ReflectedLight {
  vec3 directDiffuse;
  vec3 directSpecular;
  vec3 indirectDiffuse;
  vec3 indirectSpecular;
};

まとめると，レンダリング方程式を解くために必要なパラメータ「入射光」「幾何情報」「物体表面の材質」を渡して，「反射光」を求める処理を実装します．これは次のようになります．

// RenderEquations(RE)
void RE_Direct(const in IncidentLight directLight, const in GeometricContext geometry, const in Material material, inout ReflectedLight reflectedLight) {
  
  float dotNL = saturate(dot(geometry.normal, directLight.direction));
  vec3 irradiance = dotNL * directLight.color;
  
  // punctual light
  irradiance *= PI;
  
  reflectedLight.directDiffuse += irradiance * DiffuseBRDF(material.diffuseColor);
  reflectedLight.directSpecular += irradiance * SpecularBRDF(directLight, geometry, material.specularColor, material.specularRoughness);
}

１つ１つ見ていくと

float dotNL = saturate(dot(geometry.normal, directLight.direction));

ここでコサイン項を計算しています．

vec3 irradiance = dotNL * directLight.color;

放射輝度とコサイン項を掛けて放射照度に変換しています．

// punctual light
irradiance *= PI;

punctual light には \pi を掛けて調整します．

reflectedLight.directDiffuse += irradiance * DiffuseBRDF(material.diffuseColor);
reflectedLight.directSpecular += irradiance * SpecularBRDF(directLight, geometry, material.specularColor, material.specularRoughness);

拡散反射成分と鏡面反射成分を計算します．

そして，これを光源の数だけ処理を行います．

ReflectedLight reflectedLight = ReflectedLight(vec3(0.0), vec3(0.0), vec3(0.0), vec3(0.0));
vec3 emissive = vec3(0.0);
float opacity = 1.0;

IncidentLight directLight;

// point light
for (int i=0; i<LIGHT_MAX; ++i) {
  if (i >= numPointLights) break;
  getPointDirectLightIrradiance(pointLights[i], geometry, directLight);
  if (directLight.visible) {
    RE_Direct(directLight, geometry, material, reflectedLight);
  }
}

// spot light
for (int i=0; i<LIGHT_MAX; ++i) {
  if (i >= numSpotLights) break;
  getSpotDirectLightIrradiance(spotLights[i], geometry, directLight);
  if (directLight.visible) {
    RE_Direct(directLight, geometry, material, reflectedLight);
  }
}

// directional light
for (int i=0; i<LIGHT_MAX; ++i) {
  if (i >= numDirectionalLights) break;
  getDirectionalDirectLightIrradiance(directionalLights[i], geometry, directLight);
  RE_Direct(directLight, geometry, material, reflectedLight);
}

最後に自己放射と反射成分の総和が，最終的に出力される放射輝度になります．

vec3 outgoingLight = emissive + reflectedLight.directDiffuse + reflectedLight.directSpecular + reflectedLight.indirectDiffuse + reflectedLight.indirectSpecular;
gl_FragColor = vec4(outgoingLight, opacity);

最終コード

フラグメントシェーダの内容は以下のとおりです．

varying vec3 vViewPosition;
varying vec3 vNormal;

// uniforms
uniform float metallic;
uniform float roughness;
uniform vec3 albedo;

// defines
#define PI 3.14159265359
#define PI2 6.28318530718
#define RECIPROCAL_PI 0.31830988618
#define RECIPROCAL_PI2 0.15915494
#define LOG2 1.442695
#define EPSILON 1e-6

struct IncidentLight {
  vec3 color;
  vec3 direction;
  bool visible;
};

struct ReflectedLight {
  vec3 directDiffuse;
  vec3 directSpecular;
  vec3 indirectDiffuse;
  vec3 indirectSpecular;
};

struct GeometricContext {
  vec3 position;
  vec3 normal;
  vec3 viewDir;
};

struct Material {
  vec3 diffuseColor;
  float specularRoughness;
  vec3 specularColor;
};

// lights

bool testLightInRange(const in float lightDistance, const in float cutoffDistance) {
  return any(bvec2(cutoffDistance == 0.0, lightDistance < cutoffDistance));
}

float punctualLightIntensityToIrradianceFactor(const in float lightDistance, const in float cutoffDistance, const in float decayExponent) {
  if (decayExponent > 0.0) {
    return pow(saturate(-lightDistance / cutoffDistance + 1.0), decayExponent);
  }
  
  return 1.0;
}

struct DirectionalLight {
  vec3 direction;
  vec3 color;
};

void getDirectionalDirectLightIrradiance(const in DirectionalLight directionalLight, const in GeometricContext geometry, out IncidentLight directLight) {
  directLight.color = directionalLight.color;
  directLight.direction = directionalLight.direction;
  directLight.visible = true;
}

struct PointLight {
  vec3 position;
  vec3 color;
  float distance;
  float decay;
};

void getPointDirectLightIrradiance(const in PointLight pointLight, const in GeometricContext geometry, out IncidentLight directLight) {
  vec3 L = pointLight.position - geometry.position;
  directLight.direction = normalize(L);
  
  float lightDistance = length(L);
  if (testLightInRange(lightDistance, pointLight.distance)) {
    directLight.color = pointLight.color;
    directLight.color *= punctualLightIntensityToIrradianceFactor(lightDistance, pointLight.distance, pointLight.decay);
    directLight.visible = true;
  } else {
    directLight.color = vec3(0.0);
    directLight.visible = false;
  }
}

struct SpotLight {
  vec3 position;
  vec3 direction;
  vec3 color;
  float distance;
  float decay;
  float coneCos;
  float penumbraCos;
};

void getSpotDirectLightIrradiance(const in SpotLight spotLight, const in GeometricContext geometry, out IncidentLight directLight) {
  vec3 L = spotLight.position - geometry.position;
  directLight.direction = normalize(L);
  
  float lightDistance = length(L);
  float angleCos = dot(directLight.direction, spotLight.direction);
  
  if (all(bvec2(angleCos > spotLight.coneCos, testLightInRange(lightDistance, spotLight.distance)))) {
    float spotEffect = smoothstep(spotLight.coneCos, spotLight.penumbraCos, angleCos);
    directLight.color = spotLight.color;
    directLight.color *= spotEffect * punctualLightIntensityToIrradianceFactor(lightDistance, spotLight.distance, spotLight.decay);
    directLight.visible = true;
  } else {
    directLight.color = vec3(0.0);
    directLight.visible = false;
  }
}

// light uniforms
#define LIGHT_MAX 4
uniform DirectionalLight directionalLights[LIGHT_MAX];
uniform PointLight pointLights[LIGHT_MAX];
uniform SpotLight spotLights[LIGHT_MAX];
uniform int numDirectionalLights;
uniform int numPointLights;
uniform int numSpotLights;

// BRDFs

// Normalized Lambert
vec3 DiffuseBRDF(vec3 diffuseColor) {
  return diffuseColor / PI;
}

vec3 F_Schlick(vec3 specularColor, vec3 H, vec3 V) {
  return (specularColor + (1.0 - specularColor) * pow(1.0 - saturate(dot(V,H)), 5.0));
}

float D_GGX(float a, float dotNH) {
  float a2 = a*a;
  float dotNH2 = dotNH*dotNH;
  float d = dotNH2 * (a2 - 1.0) + 1.0;
  return a2 / (PI * d * d);
}

float G_Smith_Schlick_GGX(float a, float dotNV, float dotNL) {
  float k = a*a*0.5 + EPSILON;
  float gl = dotNL / (dotNL * (1.0 - k) + k);
  float gv = dotNV / (dotNV * (1.0 - k) + k);
  return gl*gv;
}

// Cook-Torrance
vec3 SpecularBRDF(const in IncidentLight directLight, const in GeometricContext geometry, vec3 specularColor, float roughnessFactor) {
  
  vec3 N = geometry.normal;
  vec3 V = geometry.viewDir;
  vec3 L = directLight.direction;
  
  float dotNL = saturate(dot(N,L));
  float dotNV = saturate(dot(N,V));
  vec3 H = normalize(L+V);
  float dotNH = saturate(dot(N,H));
  float dotVH = saturate(dot(V,H));
  float dotLV = saturate(dot(L,V));
  float a = roughnessFactor * roughnessFactor;

  float D = D_GGX(a, dotNH);
  float G = G_Smith_Schlick_GGX(a, dotNV, dotNL);
  vec3 F = F_Schlick(specularColor, V, H);
  return (F*(G*D))/(4.0*dotNL*dotNV+EPSILON);
}

// RenderEquations(RE)
void RE_Direct(const in IncidentLight directLight, const in GeometricContext geometry, const in Material material, inout ReflectedLight reflectedLight) {
  
  float dotNL = saturate(dot(geometry.normal, directLight.direction));
  vec3 irradiance = dotNL * directLight.color;
  
  // punctual light
  irradiance *= PI;
  
  reflectedLight.directDiffuse += irradiance * DiffuseBRDF(material.diffuseColor);
  reflectedLight.directSpecular += irradiance * SpecularBRDF(directLight, geometry, material.specularColor, material.specularRoughness);
}

void main() {
  GeometricContext geometry;
  geometry.position = -vViewPosition;
  geometry.normal = normalize(vNormal);
  geometry.viewDir = normalize(vViewPosition);
  
  Material material;
  material.diffuseColor = mix(albedo, vec3(0.0), metallic);
  material.specularColor = mix(vec3(0.04), albedo, metallic);
  material.specularRoughness = roughness;
  
  // Lighting
  
  ReflectedLight reflectedLight = ReflectedLight(vec3(0.0), vec3(0.0), vec3(0.0), vec3(0.0));
  vec3 emissive = vec3(0.0);
  float opacity = 1.0;
  
  IncidentLight directLight;
  
  // point light
  for (int i=0; i<LIGHT_MAX; ++i) {
    if (i >= numPointLights) break;
    getPointDirectLightIrradiance(pointLights[i], geometry, directLight);
    if (directLight.visible) {
      RE_Direct(directLight, geometry, material, reflectedLight);
    }
  }
  
  // spot light
  for (int i=0; i<LIGHT_MAX; ++i) {
    if (i >= numSpotLights) break;
    getSpotDirectLightIrradiance(spotLights[i], geometry, directLight);
    if (directLight.visible) {
      RE_Direct(directLight, geometry, material, reflectedLight);
    }
  }
  
  // directional light
  for (int i=0; i<LIGHT_MAX; ++i) {
    if (i >= numDirectionalLights) break;
    getDirectionalDirectLightIrradiance(directionalLights[i], geometry, directLight);
    RE_Direct(directLight, geometry, material, reflectedLight);
  }
  
  vec3 outgoingLight = emissive + reflectedLight.directDiffuse + reflectedLight.directSpecular + reflectedLight.indirectDiffuse + reflectedLight.indirectSpecular;
  
  gl_FragColor = vec4(outgoingLight, opacity);
}

📌 サンプル

動作するサンプルを用意しました．Three.js には物理ベースレンダリングを行うマテリアルがありますが，それと切り替えられるようになっています．フレネル関係で多少違うところがありますが，ほとんど変わらないはずです．

http://mebiusbox.github.io/contents/pbr/

📌 最後に

今回は拡散反射モデルと鏡面反射モデルを使った物理ベースレンダリングを実装しました．これで最低限の物理ベースレンダリングが実装されたことになったと思います．ここから，テクスチャを貼って質感を上げたり，拡散反射モデルにオーレン・ネイヤーモデルを使ったり，鏡面反射モデルは異方性を扱ったものなど多くの反射モデルがあります．また，グローバルイルミネーションでライトプローブやイラディアンスボリュームを使って間接光の拡散反射を計算した場合は reflectedLight.indirectDiffuse に，間接光の鏡面反射によって周りの景色が映り込むようなものは reflectedLight.indirectSpecular に入ってくるでしょう．シャドウマップなどで影の処理をする場合は光源から入射光を求めたときに，incidentLight.color に影の係数を掛けてから反射される輝度を求めるといったことが出来ます．

物理ベースレンダリングの基礎から実装まで，解説を書いてみましたが，いかがだったでしょうか．
これからも勉強して理解を深めつつ，記事の修正・加筆，新しい内容について書いていこうかなと思います．
少しでも参考になれば幸いです．

📌 参考文献

  • hanecci「従来のゲーム向けライト(ポイントライト, ディレクショナルライト)での BRDF の利用」＜ http://d.hatena.ne.jp/hanecci/20130604/p1 ＞（２０１７年７月２日閲覧）