こんにちわ。
Pocolです。

会社でキャラの肌が調整しずらいからどうにかしてくれ！
…と，言われてしまったので，とりあえずPre-Integrated Skinを調べてみようと思いました。

LUTテーブル作るところまでは実装してみました。
出来上がったLUTのテクスチャはこちら。

で，これを作るソースコードは以下のような感じ。

#define STBI_MSC_SECURE_CRT
#define STB_IMAGE_WRITE_IMPLEMENTATION

//-----------------------------------------------------------------------------
// Includes
//-----------------------------------------------------------------------------
#include <cstdio>
#include <stb_image_write.h>
#include <vector>
#include <asdxMath.h>

//-----------------------------------------------------------------------------
//      リニアからSRGBに変換.
//-----------------------------------------------------------------------------
inline float ToSRGB(float value)
{
    return (value <= 0.0031308) 
        ? 12.92f * value 
        : (1.0f + 0.055f) * pow(abs(value), 1.0f / 2.4f) - 0.055f;
}

//-----------------------------------------------------------------------------
//      ガウス分布.
//-----------------------------------------------------------------------------
inline float Gaussian(float v, float r)
{ return 1.0f / sqrtf(2.0f * asdx::F_PI * v) * exp(-(r * r) / (2.0f * v)); }

//-----------------------------------------------------------------------------
//      散乱計算
//-----------------------------------------------------------------------------
inline asdx::Vector3 Scatter(float r)
{
    // GPU Pro 360 Guide to Rendering, "5. Pre-Integrated Skin Shading", Appendix A.
    return 
          Gaussian(0.0064f * 1.414f, r) * asdx::Vector3(0.233f, 0.455f, 0.649f)
        + Gaussian(0.0484f * 1.414f, r) * asdx::Vector3(0.100f, 0.336f, 0.344f)
        + Gaussian(0.1870f * 1.414f, r) * asdx::Vector3(0.118f, 0.198f, 0.000f)
        + Gaussian(0.5670f * 1.414f, r) * asdx::Vector3(0.113f, 0.007f, 0.007f)
        + Gaussian(1.9900f * 1.414f, r) * asdx::Vector3(0.358f, 0.004f, 0.00001f)
        + Gaussian(7.4100f * 1.414f, r) * asdx::Vector3(0.078f, 0.00001f, 0.00001f);
}

//-------------------------------------------------------------------------------
//      Diffusionプロファイルを求める.
//-------------------------------------------------------------------------------
asdx::Vector3 IntegrateDiffuseScatterOnRing(float cosTheta, float skinRadius, int sampleCount)
{
    auto theta = acosf(cosTheta);
    asdx::Vector3 totalWeight(0.0f, 0.0f, 0.0f);
    asdx::Vector3 totalLight (0.0f, 0.0f, 0.0f);

    const auto inc = asdx::F_PI / float(sampleCount);
    auto a = -asdx::F_PIDIV2;

    while(a <= asdx::F_PIDIV2)
    {
        auto sampleAngle = theta + a;
        auto diffuse     = asdx::Saturate(cosf(sampleAngle));
        auto sampleDist  = abs(2.0f * skinRadius * sinf(a * 0.5f));
        auto weights     = Scatter(sampleDist);

        totalWeight += weights;
        totalLight  += weights * diffuse;
        a += inc;
    }

    return asdx::Vector3(
        totalLight.x / totalWeight.x,
        totalLight.y / totalWeight.y,
        totalLight.z / totalWeight.z);
}

//-----------------------------------------------------------------------------
//      ルックアップテーブル書き出し.
//-----------------------------------------------------------------------------
bool WriteLUT(const char* path, int w, int h, int s)
{
    std::vector<uint8_t> pixels;
    pixels.resize(w * h * 3);

    float stepR = 1.0f / float(h);
    float stepT = 2.0f / float(w);

    for(auto j=0; j<h; ++j)
    {
        for(auto i=0; i<w; ++i)
        {
            auto radius    = float(j + 0.5f) * stepR;
            auto curvature = 1.0f / radius;
            auto cosTheta  = -1.0f + float(i) * stepT;
            auto val       = IntegrateDiffuseScatterOnRing(cosTheta, curvature, s);

            // 書き出しを考慮して, 上下逆にして正しく出力されるようにする.
            auto idx = ((h - 1 - j) * w * 3) + (i * 3);

            pixels[idx + 0] = static_cast<uint8_t>(ToSRGB(val.x) * 255.0f + 0.5f);
            pixels[idx + 1] = static_cast<uint8_t>(ToSRGB(val.y) * 255.0f + 0.5f);
            pixels[idx + 2] = static_cast<uint8_t>(ToSRGB(val.z) * 255.0f + 0.5f);
        }
    }

    auto ret = stbi_write_tga(path, w, h, 3, pixels.data()) != 0;
    pixels.clear();

    return ret;
}

//-----------------------------------------------------------------------------
//      メインエントリーポイントです.
//-----------------------------------------------------------------------------
int main(int argc, char** argv)
{
    std::string path = "preintegrated_skin_lut.tga";
    int w = 256;
    int h = 256;
    int s = 4096;

    for(auto i=0; i<argc; ++i)
    {
        if (_stricmp(argv[i], "-w") == 0)
        {
            i++;
            auto iw = atoi(argv[i]);
            w = asdx::Clamp(iw, 1, 8192);
        }
        else if (_stricmp(argv[i], "-h") == 0)
        {
            i++;
            auto ih = atoi(argv[i]);
            h = asdx::Clamp(ih, 1, 8192);
        }
        else if (_stricmp(argv[i], "-s") == 0)
        {
            i++;
            auto is = atoi(argv[i]);
            s = asdx::Max(is, 1);
        }
    }

    return WriteLUT(path.c_str(), w, h, s) ? 0 : -1;
}

ちょっとググって調べた感じだと，https://j1jeong.wordpress.com/2018/06/20/pre-intergrated-skin-shading-in-colorspace/というBlog記事をみて，検証もせずにとりあえず鵜呑みで，sRGBで書き出してみました。
シェーダで使う場合は補正入れてリニアになるようにしてから使ってください。

良くライティング計算結果がNaNになってポストエフェクトでバグるというのが頻発していて困っていたのですが，
ようやく最近回避方法を知りました。
StackOverflowとかみていたら，

step(value, value) * value;

みたいにすると，NaNである場合には0になり，NaN以外の場合にはvalueで返すことが出来るらしい…。
と書いてあったのですが，よくよく考えてみるとstep(value, value)のところは確かにゼロとなるので問題ないなさそうに思えるのですが，NaNとの四則演算結果はNaNになるような気がします。実際に試してみたのですが，やっぱり駄目でした。

あとは，Googleのfilamentに

#define MEDIUMP_FLT_MAX    65504.0
#define saturateMediump(x) min(x, MEDIUMP_FLT_MAX)

みたいなのが書いてあったので，下記のような実装を試してみました。

static const float HALF_MAX = 65504.0;
clamp(value, 0.0f, HALF_MAX);

結果として，NaNが発生しなくなりバグが解消されるようになりました。
自分なりになぜNaNが発生しなくなるのかを考えてみたのですが，NaNの性質として !=以外の演算は常にfalse, !=演算は常にtrueという性質があるようです。よって，clamp()メソッドを使うと比較演算が実行されるので，この結果としてNaNをはじいて指定範囲内に収まると考えれば非常に納得がいきます。ここではclamp()を試しましたが，同様の理屈が通るのであればmax()などでも代用できるはずです。

さて，特に物理ベースレンダリングに移行する際に，GGXモデルなどを使うことが多いと思うのですが，定式化されているものが除算を含む形で定義されているので，実装する際に除算をしなくてはならないということが発生します。
GGXの計算する際に，分母がゼロになることにより，ゼロ除算が発生し，計算結果がNaNになる可能性があります。
また，pow()関数を使用している場合は，MSDNのドキュメント にも書いてあるように第1引数が0未満の場合はNaNになり，第1引数・第2引数がともにゼロの場合はハードウェア依存によりNaNが発生する可能性があります。

昔からよくある方法としてゼロ除算を発生させないために，例えば1e-7fなどのような小さな値を分母側に足しておき，ゼロにならないようにオフセットを掛けるという回避方法もあるのですが，この方法は除算には適用できるのだけれども，pow()関数には適用できません。

GGXを用いたライティング計算の場合は，ゼロ除算が発生する可能性があるため最初はこのオフセットを足し込むというやり方を取っていたのですが，これをやってしまうと計算結果が変わってしまい，ハイライトが鈍くなるなど見た目上に影響を及ぼすことが分かりました。
PBRやっているはずなのに，何故かなまった感じの画が出る場合には，下駄をはかせてハイライトを潰してしまっていないかどうか確認するのを強くお勧めします。
下駄をはかせている箇所を下駄を取って普通に除算し，最後に上記のようなNaNをはじく処理をいれることで，ハイライトがなまったり，ライティング以降のレンダリングパスでバグるという問題を解決することが出来ます。また，pow()関数を用いる場合にも適用でき，NaNをなくすることができるので，NaNで困っている方はぜひ試してみてください。

※もっとより良い方法ご存知であれば，是非教えてください！