图:大理石参数渲染结果
概述
我在刚开始接触十多个参数的Blender Disney BSDF时就非常麻爪,很多讲解的文章都停留在单个参数讲解的层面上,也搜索不到现成的参数参考。然后我看到了毛星云(RIP🙏🏼)在PBR白皮书的Disney Principled BSDF一节中,提到迪士尼根据对MERL100数据库的观察与理论分析提出了该模型。
有了MERL100数据库,就可以计算获得更接近真实的材质参数。换句话说,从MERL100计算Disney Principled BSDF的参数也可以称得上是原汤化原食了。
本文中,Disney Principled BSDF的代码来自github的源码,参数与blender内置的Disney Principled BRDF有一定差异,主要表现为Blender细化了subsurface颜色、IOR和各向异性细化为过滤和旋转,对材质有些微影响。
MERL100的数据由Wojciech Matusik 、 Hanspeter Pfister 、 Matt Brand 和 Leonard McMillan测得,数据文件来自知乎@xkl的分享。数据文件由binary数据和c++编译器组成。进行c++编译之后可以得到RGB值。
本文计算依据的公式是BSDF * Light(光照参数) = MERL,在观察结果时发现MERL颜色偏黑,就加入了MERL_fix参数提亮画面,最终公式为BSDF * Light = MERL * MERL_fix。
渲染结果如标题图所示,但是相同材质在不同环境光的表现可能大相径庭,MERL100算出来的参数只能说是一根比较好用的拐杖吧。
零、MERL100数据库分析
1.MERL100材质库包含的材质
beige-fabric | 米色织物 | 织物 | pink-plastic | 粉色塑料 | 塑料 |
black-fabric | 黑色织物 | 织物 | polyethylene | 聚乙烯 | 塑料 |
blue-fabric | 蓝色织物 | 织物 | polyurethane-foam | 聚氨酯泡沫 | 塑料 |
dark-specular-fabric | 高反射织物 | 织物 | pure-rubber | 天然橡胶 | 塑料 |
green-fabric | 绿色织物 | 织物 | pvc | 聚氯乙烯 | 塑料 |
light-brown-fabric | 浅棕色织物 | 织物 | red-phenolic | 红色酚醛塑料 | 塑料 |
nylon | 尼龙 | 织物 | red-plastic | 红色塑料 | 塑料 |
pink-fabric | 粉色织物 | 织物 | red-specular-plastic | 镜面红色塑料 | 塑料 |
pink-fabric2 | 粉色织物 | 织物 | specular-black-phenolic | 镜面黑色酚醛塑料 | 塑料 |
pink-felt | 粉色毛毡 | 织物 | specular-blue-phenolic | 镜面蓝色酚醛塑料 | 塑料 |
red-fabric | 红色织物 | 织物 | specular-green-phenolic | 镜面绿色酚醛塑料 | 塑料 |
red-fabric2 | 红色织物 | 织物 | specular-maroon-phenolic | 镜面栗色酚醛塑料 | 塑料 |
teflon | 聚四氟乙烯 | 织物 | specular-orange-phenolic | 镜面橘色酚醛塑料 | 塑料 |
white-fabric | 白色织物 | 织物 | specular-red-phenolic | 镜面红色酚醛塑料 | 塑料 |
white-fabric2 | 白色织物 | 织物 | specular-violet-phenolic | 镜面紫罗兰色酚醛塑料 | 塑料 |
blue-acrylic | 蓝丙烯 | 涂料 | specular-white-phenolic | 镜面白色酚醛塑料 | 塑料 |
blue-metallic-paint | 蓝色金属漆 | 涂料 | specular-yellow-phenolic | 镜面黄色酚醛塑料 | 塑料 |
blue-metallic-paint2 | 蓝色金属漆2 | 涂料 | violet-acrylic | 紫色丙烯酸 | 塑料 |
color-changing-paint1 | 变色颜料1 | 涂料 | violet-rubber | 紫色橡胶 | 塑料 |
color-changing-paint2 | 变色颜料2 | 涂料 | white-diffuse-bball | 白色漫反射球 | 塑料 |
color-changing-paint3 | 变色颜料3 | 涂料 | yellow-matte-plastic | 黄色哑光塑料 | 塑料 |
dark-blue-paint | 深蓝涂料 | 涂料 | yellow-phenolic | 黄色酚醛塑料 | 塑料 |
dark-red-paint | 深红涂料 | 涂料 | yellow-plastic | 黄色塑料 | 塑料 |
gold-metallic-paint | 金色金属漆 | 涂料 | cherry-235 | cherry235(木材清漆 | 木材 |
gold-metallic-paint2 | 金色金属漆2 | 涂料 | colonial-maple-223 | CM223(木材清漆 | 木材 |
gold-metallic-paint3 | 金色金属漆3 | 涂料 | fruitwood-241 | fw241(木材清漆 | 木材 |
gold-paint | 金色涂料 | 涂料 | ipswich-pine-221 | ip221(木材清漆 | 木材 |
green-latex | 绿色乳胶 | 涂料 | natural-209 | n209(木材清漆 | 木材 |
green-metallic-paint | 绿色金属涂料 | 涂料 | pickled-oak-260 | po260(木材清漆 | 木材 |
green-metallic-paint2 | 绿色金属涂料2 | 涂料 | special-walnut-224 | sw224(木材清漆 | 木材 |
light-red-paint | 浅红色涂料 | 涂料 | aventurnine | 东陵玉 | 矿物 |
orange-paint | 橘色涂料 | 涂料 | black-obsidian | 黑曜石 | 矿物 |
pearl-paint | 珠光漆 | 涂料 | pink-jasper | 粉碧玉 | 矿物 |
purple-paint | 紫色涂料 | 涂料 | silicon-nitrade | 氮化硅 | 矿物 |
red-metallic-paint | 红色金属漆 | 涂料 | white-marble | 白色大理石 | 矿物 |
silver-metallic-paint | 银色金属漆 | 涂料 | alum-bronze | 铝青铜 | 金属 |
silver-metallic-paint2 | 银色金属漆 | 涂料 | alumina-oxide | 氧化铝 | 金属 |
silver-paint | 银色涂料 | 涂料 | aluminium | 铝 | 金属 |
white-acrylic | 白色丙烯酸 | 涂料 | black-oxidized-steel | 染黑铁 | 金属 |
white-paint | 白色涂料 | 涂料 | brass | 黄铜 | 金属 |
yellow-paint | 黄色涂料 | 涂料 | chrome-steel | 铬钢 | 金属 |
black-phenolic | 黑色酚醛塑料 | 塑料 | chrome | 铬 | 金属 |
black-soft-plastic | 黑色柔性塑料 | 塑料 | grease-covered-steel | 涂油钢 | 金属 |
blue-rubber | 蓝色橡胶 | 塑料 | hematite | 赤铁矿 | 金属 |
delrin | 缩醛树脂 | 塑料 | nickel | 镍 | 金属 |
gray-plastic | 灰色塑料 | 塑料 | ss440 | ss440碳素钢 | 金属 |
green-acrylic | 绿色丙烯酸 | 塑料 | steel | 钢 | 金属 |
green-plastic | 绿色塑料 | 塑料 | tungsten-carbide | 碳化钨 | 金属 |
maroon-plastic | 栗色塑料 | 塑料 | two-layer-gold | 两层金 | 金属 |
neoprene-rubber | 氯丁橡胶 | 塑料 | two-layer-silver | 两层银 | 金属 |
2.cpp编译结果和改写
BRDF-pic-Read.cpp运行结果为RGB值文件,符合Disney Principled BxDF所需的结果,Disney Principled BxDF的输入参数NLV在cpp计算过程中出现过,可通过改写cpp得到:
// Copyright 2005 Mitsubishi Electric Research Laboratories All Rights Reserved.
// Permission to use, copy and modify this software and its documentation without
// fee for educational, research and non-profit purposes, is hereby granted, provided
// that the above copyright notice and the following three paragraphs appear in all copies.
// To request permission to incorporate this software into commercial products contact:
// Vice President of Marketing and Business Development;
// Mitsubishi Electric Research Laboratories (MERL), 201 Broadway, Cambridge, MA 02139 or
// <license@merl.com>.
// IN NO EVENT SHALL MERL BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL,
// OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
// ITS DOCUMENTATION, EVEN IF MERL HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
// MERL SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED
// HEREUNDER IS ON AN "AS IS" BASIS, AND MERL HAS NO OBLIGATIONS TO PROVIDE MAINTENANCE, SUPPORT,
// UPDATES, ENHANCEMENTS OR MODIFICATIONS.
#include "stdlib.h";
#include "math.h";
#include <cstdio>;
#define BRDF_SAMPLING_RES_THETA_H 90
#define BRDF_SAMPLING_RES_THETA_D 90
#define BRDF_SAMPLING_RES_PHI_D 360
#define RED_SCALE (1.0/1500.0)
#define GREEN_SCALE (1.15/1500.0)
#define BLUE_SCALE (1.66/1500.0)
#define M_PI 3.1415926535897932384626433832795
// cross product of two vectors
void cross_product (double* v1, double* v2, double* out)
{
out[0] = v1[1]*v2[2] - v1[2]*v2[1];
out[1] = v1[2]*v2[0] - v1[0]*v2[2];
out[2] = v1[0]*v2[1] - v1[1]*v2[0];
}
// normalize vector
void normalize(double* v)
{
// normalize
double len = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
v[0] = v[0] / len;
v[1] = v[1] / len;
v[2] = v[2] / len;
}
// rotate vector along one axis
void rotate_vector(double* vector, double* axis, double angle, double* out)
{
double temp;
double cross[3];
double cos_ang = cos(angle);
double sin_ang = sin(angle);
out[0] = vector[0] * cos_ang;
out[1] = vector[1] * cos_ang;
out[2] = vector[2] * cos_ang;
temp = axis[0]*vector[0]+axis[1]*vector[1]+axis[2]*vector[2];
temp = temp*(1.0-cos_ang);
out[0] += axis[0] * temp;
out[1] += axis[1] * temp;
out[2] += axis[2] * temp;
cross_product (axis,vector,cross);
out[0] += cross[0] * sin_ang;
out[1] += cross[1] * sin_ang;
out[2] += cross[2] * sin_ang;
}
// convert standard coordinates to half vector/difference vector coordinates
void std_coords_to_half_diff_coords(double theta_in, double fi_in, double theta_out, double fi_out,
double& theta_half,double& fi_half,double& theta_diff,double& fi_diff )
{
// compute in vector
double in_vec_z = cos(theta_in);
double proj_in_vec = sin(theta_in);
double in_vec_x = proj_in_vec*cos(fi_in);
double in_vec_y = proj_in_vec*sin(fi_in);
double in[3]= {in_vec_x,in_vec_y,in_vec_z};
normalize(in);
// compute out vector
double out_vec_z = cos(theta_out);
double proj_out_vec = sin(theta_out);
double out_vec_x = proj_out_vec*cos(fi_out);
double out_vec_y = proj_out_vec*sin(fi_out);
double out[3]= {out_vec_x,out_vec_y,out_vec_z};
normalize(out);
// compute halfway vector
double half_x = (in_vec_x + out_vec_x)/2.0f;
double half_y = (in_vec_y + out_vec_y)/2.0f;
double half_z = (in_vec_z + out_vec_z)/2.0f;
double half[3] = {half_x,half_y,half_z};
normalize(half);
printf("[%f %f %f], [%f %f %f], [%f %f %f]\n", (float)in[0], (float)in[1], (float)in[2], (double)out[0], (double)out[1], (double)out[2], (double)half[0], (double)half[1], (double)half[2]);
// compute theta_half, fi_half
theta_half = acos(half[2]);
fi_half = atan2(half[1], half[0]);
double bi_normal[3] = {0.0, 1.0, 0.0};
double normal[3] = { 0.0, 0.0, 1.0 };
double temp[3];
double diff[3];
// compute diff vector
rotate_vector(in, normal , -fi_half, temp);
rotate_vector(temp, bi_normal, -theta_half, diff);
// compute theta_diff, fi_diff
theta_diff = acos(diff[2]);
fi_diff = atan2(diff[1], diff[0]);
}
// Lookup theta_half index
// This is a non-linear mapping!
// In: [0 .. pi/2]
// Out: [0 .. 89]
inline int theta_half_index(double theta_half)
{
if (theta_half <= 0.0)
return 0;
double theta_half_deg = ((theta_half / (M_PI/2.0))*BRDF_SAMPLING_RES_THETA_H);
double temp = theta_half_deg*BRDF_SAMPLING_RES_THETA_H;
temp = sqrt(temp);
int ret_val = (int)temp;
if (ret_val < 0) ret_val = 0;
if (ret_val >= BRDF_SAMPLING_RES_THETA_H)
ret_val = BRDF_SAMPLING_RES_THETA_H-1;
return ret_val;
}
// Lookup theta_diff index
// In: [0 .. pi/2]
// Out: [0 .. 89]
inline int theta_diff_index(double theta_diff)
{
int tmp = int(theta_diff / (M_PI * 0.5) * BRDF_SAMPLING_RES_THETA_D);
if (tmp < 0)
return 0;
else if (tmp < BRDF_SAMPLING_RES_THETA_D - 1)
return tmp;
else
return BRDF_SAMPLING_RES_THETA_D - 1;
}
// Lookup phi_diff index
inline int phi_diff_index(double phi_diff)
{
// Because of reciprocity, the BRDF is unchanged under
// phi_diff -> phi_diff + M_PI
if (phi_diff < 0.0)
phi_diff += M_PI;
// In: phi_diff in [0 .. pi]
// Out: tmp in [0 .. 179]
int tmp = int(phi_diff / M_PI * BRDF_SAMPLING_RES_PHI_D / 2);
if (tmp < 0)
return 0;
else if (tmp < BRDF_SAMPLING_RES_PHI_D / 2 - 1)
return tmp;
else
return BRDF_SAMPLING_RES_PHI_D / 2 - 1;
}
// Given a pair of incoming/outgoing angles, look up the BRDF.
void lookup_brdf_val(double* brdf, double theta_in, double fi_in,
double theta_out, double fi_out,
double& red_val,double& green_val,double& blue_val)
{
// Convert to halfangle / difference angle coordinates
double theta_half, fi_half, theta_diff, fi_diff;
std_coords_to_half_diff_coords(theta_in, fi_in, theta_out, fi_out,
theta_half, fi_half, theta_diff, fi_diff);
// Find index.
// Note that phi_half is ignored, since isotropic BRDFs are assumed
int ind = phi_diff_index(fi_diff) +
theta_diff_index(theta_diff) * BRDF_SAMPLING_RES_PHI_D / 2 +
theta_half_index(theta_half) * BRDF_SAMPLING_RES_PHI_D / 2 *
BRDF_SAMPLING_RES_THETA_D;
//printf("%f\n", (float)ind);
red_val = brdf[ind] * RED_SCALE;
green_val = brdf[ind + BRDF_SAMPLING_RES_THETA_H*BRDF_SAMPLING_RES_THETA_D*BRDF_SAMPLING_RES_PHI_D/2] * GREEN_SCALE;
blue_val = brdf[ind + BRDF_SAMPLING_RES_THETA_H*BRDF_SAMPLING_RES_THETA_D*BRDF_SAMPLING_RES_PHI_D] * BLUE_SCALE;
//printf("%f %f %f %f\n", (float)ind, (float)brdf[ind], (float)RED_SCALE, (float)red_val);
if (red_val < 0.0 || green_val < 0.0 || blue_val < 0.0)
fprintf(stderr, "Below horizon.\n");
}
// Read BRDF data
bool read_brdf(const char *filename, double* &brdf)
{
FILE *f = fopen(filename, "rb");
if (!f)
return false;
int dims[3];
fread(dims, sizeof(int), 3, f);
int n = dims[0] * dims[1] * dims[2];
if (n != BRDF_SAMPLING_RES_THETA_H *
BRDF_SAMPLING_RES_THETA_D *
BRDF_SAMPLING_RES_PHI_D / 2)
{
fprintf(stderr, "Dimensions don't match\n");
fclose(f);
return false;
}
brdf = (double*) malloc (sizeof(double)*3*n);
fread(brdf, sizeof(double), 3*n, f);
fclose(f);
return true;
}
int main(int argc, char *argv[])
{
const char *filename = argv[1];
double* brdf;
// read brdf
if (!read_brdf(filename, brdf))
{
fprintf(stderr, "Error reading %s\n", filename);
exit(1);
}
// print out a 16x64x16x64 table table of BRDF values
const int n = 16;
for (int i = 0; i < n; i++)
{
double theta_in = i * 0.5 * M_PI / n;
for (int j = 0; j < 4*n; j++)
{
double phi_in = j * 2.0 * M_PI / (4*n);
for (int k = 0; k < n; k++)
{
double theta_out = k * 0.5 * M_PI / n;
for (int l = 0; l < 4*n; l++)
{
double phi_out = l * 2.0 * M_PI / (4*n);
double red,green,blue;
lookup_brdf_val(brdf, theta_in, phi_in, theta_out, phi_out, red, green, blue);
printf("%f %f %f\n", (float)red, (float)green, (float)blue);
}
}
}
}
return 0;
}
3.MERL_Fix
BRDF-pic-Read.cpp编译出的rgb文件,由PIL转为图片后发现颜色过暗,引入Merl_Fix使计算结果更符合人眼观察。(例如此材料的MREL_Fix为2000)
#-*- coding:utf-8 -*-
from PIL import Image
import re
x = 1024 #x坐标 通过对txt里的行数进行整数分解
y = 1024 #y坐标 x*y = 行数
o = 2000 #MERL_Fix参数
im = Image.new("RGB",(x,y))#创建图片
file = open('fruitwood-241-color.txt') #打开rbg值文件
#通过一个个rgb点生成图片
for i in range(0,x):
for j in range(0,y):
line = file.readline().replace('[','').replace(']','')#获取一行
rgb = line.split(",")#分离rgb
#print(type(rgb[0]))
im.putpixel((i,j),(int(float(rgb[0])*o),int(float(rgb[1])*o),int(float(rgb[2])*o)))#rgb转化为像素
im.show()
一、估算Light系数
通过阅读论文,可以得到可见光发射和灯距:
通过积分估算灯亮度为20W,灯距50cm,依据条件进行blender建模:
将设定的材质代入Disney Principled BSDF进行计算,与blender渲染的RGB相除,可得light系数为600.
二、计算参数
1.Sympy直接解方程
一开始的思路肯定是怎么省事怎么来,想着直接化简到Sympy里解,然后对Disney BSDF进行化简、并把RGB三维拆分为R+G+B三个一维就得到了
ax_ = max(.001, (roughness)**2/math.sqrt(1-anistrophic*.9))
ay_ = max(.001, (roughness)**2*math.sqrt(1-anistrophic*.9))
Cdlum_ = 0.3*baseColor0**2.2 + 0.6*baseColor1**2.2 + 0.1*baseColor2**2.2
BSDF0 = (ax**3*ay**3/(PI*(NdotL + sqrt(LdotX**2*ax**2 + LdotY**2*ay**2 + NdotL**2))*(NdotV + sqrt(NdotV**2 + VdotX**2*ax**2 + VdotY**2*ay**2))*(HdotX**2*ay**2 + HdotY**2*ax**2 + NdotH**2*ax**2*ay**2)**2))*((specular*.08*((1-specularTint)+baseColor0**2.2/Cdlum*specularTint)*(1-metalic) + baseColor0**2.2*metalic)*(1-FH)+FH) + .25*clearCoat*(((.1*(1-clearCoatGloss)+.001*clearCoatGloss)**2-1) / (PI*sympy.log((.1*(1-clearCoatGloss)+.001*clearCoatGloss)**2)/sympy.log(10)*(1 + ((.1*(1-clearCoatGloss)+.001*clearCoatGloss)**2-1)*NdotH*NdotH)))*( .04*(1-FH) + 1*FH)*(1 / (NdotL + sqrt(0.25*0.25 + NdotL*NdotL - 0.25*0.25*NdotL*NdotL))* 1 / (NdotV + sqrt(0.25*0.25 + NdotV*NdotV - 0.25*0.25*NdotV*NdotV)))+ (1-metalic) * (( (1/PI) * (((1-FL)+(0.5 + 2 * LdotH*LdotH * roughness) * FL) * ((1-FV)+(0.5 + 2 * LdotH*LdotH * roughness) * FV) * (1-subsurface)+ (1.25 * ((1-FL+LdotH*LdotH*roughness*FL) * (1-FV+LdotH*LdotH*roughness*FV) * (1 / (NdotL + NdotV) - .5) + .5)) * subsurface)*baseColor0**2.2)+(FH * sheen * ((1-sheenTint)+baseColor0**2.2/Cdlum * sheenTint)))
模型参数如下:
color baseColor
float metallic
float subsurface
float specular
float roughness
float specularTint
float anisotropic
float sheen
float sheenTint
float clearcoat
float clearcoatGloss
作为自变量的方向向量:
def SchlickFresnel(u):
yy = [1-u, 0, 1]
yy.sort()
m = yy[1]
# m = clamp(1-u, 0, 1);
m2 = m*m
return m2*m2*m # pow(m,5)
def normalize(x):
return x/np.linalg.norm(x)
L = normalize(np.array([]))
V = normalize(np.array([]))
N = normalize(np.array([]))
#H = normalize(L+V)
X = np.array([1,0,0])
Y = np.array([0,1,0])
Z = np.array([0,0,1])
NdotL = np.dot(N,L)
NdotV = np.dot(N,V)
FL = SchlickFresnel(NdotL)
FV = SchlickFresnel(NdotV)
H_origion = L+V
H = H_origion/np.linalg.norm(H_origion)
NdotH = np.dot(N,H);
LdotH = np.dot(L,H);
FH = SchlickFresnel(LdotH)
VdotX = np.dot(V,X)
VdotY = np.dot(V,Y)
HdotX = np.dot(H,X)
HdotY = np.dot(H,Y)
LdotX = np.dot(L,X)
LdotY = np.dot(L,Y)
看着很复杂,果然个把小时也没算出个结果,又申请不到超级计算机,还是得换个方法。
2.随机数解方程
比较简单粗暴,直接用随机数去碰一个误差可以接受的解。然而,算出来的参数五颜六色,字面意思上的五颜六色,数值上误差很小,但视觉上根本就不是同一个东西。
baseColor0 = np.random.uniform(0,1,1000000)
baseColor1 = np.random.uniform(0,1,1000000)
baseColor2 = np.random.uniform(0,1,1000000)
metalic = np.random.uniform(0,1,1000000)
subsurface = np.random.uniform(0,1,1000000)
specular = np.random.uniform(0,1,1000000)
specularTint = np.random.uniform(0,1,1000000)
roughness = np.random.uniform(0,1,1000000)
anistrophic = np.random.uniform(0,1,1000000)
sheen = np.random.uniform(0,1,1000000)
sheenTint = np.random.uniform(0,1,1000000)
clearCoat = np.random.uniform(0,1,1000000)
clearCoatGloss = np.random.uniform(0,1,1000000)
3.优化随机数,增加随机条件
MERL100采用白光,不改变材质的颜色表现, 我根据HSV颜色锁定H色相,减少随机过程的自由度,使得计算结果更加实用。
def HSV2RGB(H,S,V):
global R,G,B
C = V * S
X = C * (1- abs(((H/60)%2) -1 ))
m = V - C
if H < 60:
R_ = C
G_ = X
B_ = 0
elif (H>=60) and (H<120):
R_ = X
G_ = C
B_ = 0
elif (H>=120) and (H<180):
R_ = 0
G_ = C
B_ = X
elif (H>=180) and (H<240):
R_ = 0
G_ = X
B_ = C
elif (H>=240) and (H<300):
R_ = X
G_ = 0
B_ = C
else:
R_ = C
G_ = 0
B_ = X
R = (R_ + m )
G = (G_ + m )
B = (B_ + m )
return R,G,B
H = 29
S = np.random.uniform(0,1,1000000)
V = np.random.uniform(0,1,1000000)
metalic = np.random.uniform(0,1,1000000)
subsurface = np.random.uniform(0,1,1000000)
specular = np.random.uniform(0,1,1000000)
specularTint = np.random.uniform(0,1,1000000)
roughness = np.random.uniform(0,1,1000000)
anistrophic = np.random.uniform(0,1,1000000)
sheen = np.random.uniform(0,1,1000000)
sheenTint = np.random.uniform(0,1,1000000)
clearCoat = np.random.uniform(0,1,1000000)
clearCoatGloss = np.random.uniform(0,1,1000000)
三、最终代码
#统计行数-均匀间隔选15个
import codecs
import numpy as np
import math
def HSV2RGB(H,S,V):
global R,G,B
C = V * S
X = C * (1- abs(((H/60)%2) -1 ))
m = V - C
if H < 60:
R_ = C
G_ = X
B_ = 0
elif (H>=60) and (H<120):
R_ = X
G_ = C
B_ = 0
elif (H>=120) and (H<180):
R_ = 0
G_ = C
B_ = X
elif (H>=180) and (H<240):
R_ = 0
G_ = X
B_ = C
elif (H>=240) and (H<300):
R_ = X
G_ = 0
B_ = C
else:
R_ = C
G_ = 0
B_ = X
R = (R_ + m )
G = (G_ + m )
B = (B_ + m )
return R,G,B
result = []
file = open ("fruitwood-241-all.txt", "r")
lines = len(file.readlines())
#print(lines)
n0 = 125
if int((lines - n0) / 30)%2 == 0:
gap = int((lines - n0) / 30)
else:
gap = int((lines - n0) / 30 - 1)
file.close()
#print(gap)
lines_result_ = []
for x in range(15):
with codecs.open('./fruitwood-241-all.txt', 'r', 'gb18030') as infile:
m = infile.readlines()[int(n0) + int(x*gap)]
#print((m))
lines_result_.append(m)
#lines_cal.append(m[0])
#print(lines_result_)
lines_result = []
for r in lines_result_:
r_ = r.replace('\r\n','').replace(' ',',')
r__ = r_.split(',')
r__ = np.float_(r__)
rNew = np.array(r__)
lines_result.append((rNew))
#print(lines_result)
lines_cal_ = []
for x in range(15):
with codecs.open('./fruitwood-241-all.txt', 'r', 'gb18030') as infile:
n = infile.readlines()[int(n0) + int(x*gap) -1]
#print((n))
lines_cal_.append(n)
#print(lines_cal_)
lines_cal = []
for f in lines_cal_:
f__ = f.split(',')
#print(len(f__))
for z in f__[0:3]:
#print(type(z))
z_ = z.replace('[','').replace(' ', ',').replace(']','').replace('\r\n','')
z__ = z_.split(',')
if(len(z__) == 4):
z__ = z__[1::]
# print(z__)
z__ = np.float_(z__)
f__.append(z__)
fnew = f__[3::]
#print(fnew)
lines_cal.append((fnew))
#print(lines_cal)
from sympy import *
import sympy
def SchlickFresnel(u):
yy = [1-u, 0, 1]
yy.sort()
m = yy[1]
# m = clamp(1-u, 0, 1);
m2 = m*m
return m2*m2*m # pow(m,5)
#算法改写
#math.sqrt→sqrt
#log10()→log()/log10
#baseColor[0]→baseColor0
#需求解的参数
roughness = symbols('roughness')
anistrophic = symbols('anistrophic')
specular = symbols('specular')
specularTint = symbols('specularTint')
baseColor0 = symbols('baseColor0')
baseColor1 = symbols('baseColor1')
baseColor2 = symbols('baseColor2')
metalic = symbols('metalic')
subsurface = symbols('subsurface')
sheen = symbols('sheen')
sheenTint = symbols('sheenTint')
clearCoat = symbols('clearCoat')
clearCoatGloss = symbols('clearCoatGloss')
ax = symbols('ax')#max(.001, (roughness)**2/math.sqrt(1-anistrophic*.9))
ay = symbols('ay')#max(.001, (roughness)**2*math.sqrt(1-anistrophic*.9))
Cdlum = symbols('Cdlum')#(0.3*Cdlin[0] + .6*Cdlin[1] + .1*Cdlin[2])>0?Cdlin/Cdlum:np.array([1,1,1])
#需替换 .subs()
PI = symbols('PI')
NdotL = symbols('NdotL')
LdotX = symbols('LdotX')
LdotY = symbols('LdotY')
NdotV = symbols('NdotV')
VdotX = symbols('VdotX')
VdotY = symbols('VdotY')
HdotX = symbols('HdotX')
HdotY = symbols('HdotY')
NdotH = symbols('NdotH')
LdotH = symbols('LdotH')
FL = symbols('FL')
FV = symbols('FV')
FH = symbols('FH')
BSDF0 = (ax**3*ay**3/(PI*(NdotL + sqrt(LdotX**2*ax**2 + LdotY**2*ay**2 + NdotL**2))*(NdotV + sqrt(NdotV**2 + VdotX**2*ax**2 + VdotY**2*ay**2))*(HdotX**2*ay**2 + HdotY**2*ax**2 + NdotH**2*ax**2*ay**2)**2))*((specular*.08*((1-specularTint)+baseColor0**2.2/Cdlum*specularTint)*(1-metalic) + baseColor0**2.2*metalic)*(1-FH)+FH) + .25*clearCoat*(((.1*(1-clearCoatGloss)+.001*clearCoatGloss)**2-1) / (PI*sympy.log((.1*(1-clearCoatGloss)+.001*clearCoatGloss)**2)/sympy.log(10)*(1 + ((.1*(1-clearCoatGloss)+.001*clearCoatGloss)**2-1)*NdotH*NdotH)))*( .04*(1-FH) + 1*FH)*(1 / (NdotL + sqrt(0.25*0.25 + NdotL*NdotL - 0.25*0.25*NdotL*NdotL))* 1 / (NdotV + sqrt(0.25*0.25 + NdotV*NdotV - 0.25*0.25*NdotV*NdotV)))+ (1-metalic) * (( (1/PI) * (((1-FL)+(0.5 + 2 * LdotH*LdotH * roughness) * FL) * ((1-FV)+(0.5 + 2 * LdotH*LdotH * roughness) * FV) * (1-subsurface)+ (1.25 * ((1-FL+LdotH*LdotH*roughness*FL) * (1-FV+LdotH*LdotH*roughness*FV) * (1 / (NdotL + NdotV) - .5) + .5)) * subsurface)*baseColor0**2.2)+(FH * sheen * ((1-sheenTint)+baseColor0**2.2/Cdlum * sheenTint)))
equation = []
solution = []
merl_fix = 2000
light_fix = 600
for x in range(14):
cals = lines_cal[x]
N = cals[2]
H = cals[2]
L = cals[0]
V = cals[1]
X = np.array([1,0,0])
Y = np.array([0,1,0])
Z = np.array([0,0,1])
PI_ = 3.14159265358979323846
NdotL_ = np.dot(N,L)
LdotX_ = np.dot(L,X)
LdotY_ = np.dot(L,Y)
NdotV_ = np.dot(N,V)
VdotX_ = np.dot(V,X)
VdotY_ = np.dot(V,Y)
HdotX_ = np.dot(H,X)
HdotY_ = np.dot(H,Y)
NdotH_ = np.dot(N,H)
LdotH_ = np.dot(L,H)
FL_ = SchlickFresnel(NdotL_)
FV_ = SchlickFresnel(NdotV_)
FH_ = SchlickFresnel(LdotH_)
BSDF0 = BSDF0.subs(PI,PI_).subs(NdotL,NdotL_).subs(LdotX,LdotX_).subs(LdotY,LdotY_).subs(NdotV,NdotV_).subs(VdotX,VdotX_).subs(VdotY,VdotY_).subs(HdotX,HdotX_).subs(HdotY,HdotY_).subs(NdotH,NdotH_).subs(LdotH,LdotH_).subs(FL,FL_).subs(FV,FV_).subs(FH,FH_)
BSDF0 = BSDF0 - lines_result[x][0]*merl_fix/light_fix/255
equation.append(BSDF0)
#print((FV_))
#print((equation[1]))
f0 = lambdify([baseColor0,baseColor1,baseColor2,metalic,subsurface,specular,specularTint, roughness,anistrophic,sheen,sheenTint,clearCoat,clearCoatGloss,ax,ay,Cdlum],equation[0],'numpy')
f1 = lambdify([baseColor0,baseColor1,baseColor2,metalic,subsurface,specular,specularTint, roughness,anistrophic,sheen,sheenTint,clearCoat,clearCoatGloss,ax,ay,Cdlum],equation[1],'numpy')
f2 = lambdify([baseColor0,baseColor1,baseColor2,metalic,subsurface,specular,specularTint, roughness,anistrophic,sheen,sheenTint,clearCoat,clearCoatGloss,ax,ay,Cdlum],equation[2],'numpy')
f3 = lambdify([baseColor0,baseColor1,baseColor2,metalic,subsurface,specular,specularTint, roughness,anistrophic,sheen,sheenTint,clearCoat,clearCoatGloss,ax,ay,Cdlum],equation[3],'numpy')
f4 = lambdify([baseColor0,baseColor1,baseColor2,metalic,subsurface,specular,specularTint, roughness,anistrophic,sheen,sheenTint,clearCoat,clearCoatGloss,ax,ay,Cdlum],equation[4],'numpy')
f5 = lambdify([baseColor0,baseColor1,baseColor2,metalic,subsurface,specular,specularTint, roughness,anistrophic,sheen,sheenTint,clearCoat,clearCoatGloss,ax,ay,Cdlum],equation[5],'numpy')
f6 = lambdify([baseColor0,baseColor1,baseColor2,metalic,subsurface,specular,specularTint, roughness,anistrophic,sheen,sheenTint,clearCoat,clearCoatGloss,ax,ay,Cdlum],equation[6],'numpy')
f7 = lambdify([baseColor0,baseColor1,baseColor2,metalic,subsurface,specular,specularTint, roughness,anistrophic,sheen,sheenTint,clearCoat,clearCoatGloss,ax,ay,Cdlum],equation[7],'numpy')
f8 = lambdify([baseColor0,baseColor1,baseColor2,metalic,subsurface,specular,specularTint, roughness,anistrophic,sheen,sheenTint,clearCoat,clearCoatGloss,ax,ay,Cdlum],equation[8],'numpy')
f9 = lambdify([baseColor0,baseColor1,baseColor2,metalic,subsurface,specular,specularTint, roughness,anistrophic,sheen,sheenTint,clearCoat,clearCoatGloss,ax,ay,Cdlum],equation[9],'numpy')
f10 = lambdify([baseColor0,baseColor1,baseColor2,metalic,subsurface,specular,specularTint, roughness,anistrophic,sheen,sheenTint,clearCoat,clearCoatGloss,ax,ay,Cdlum],equation[10],'numpy')
f11 = lambdify([baseColor0,baseColor1,baseColor2,metalic,subsurface,specular,specularTint, roughness,anistrophic,sheen,sheenTint,clearCoat,clearCoatGloss,ax,ay,Cdlum],equation[11],'numpy')
f13 = lambdify([baseColor0,baseColor1,baseColor2,metalic,subsurface,specular,specularTint, roughness,anistrophic,sheen,sheenTint,clearCoat,clearCoatGloss,ax,ay,Cdlum],equation[13],'numpy')
#yyy = f(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
#print(yyy)
#yyy = o
import numpy as np
import math
def normalize(x):
return x/np.linalg.norm(x)
def SchlickFresnel(u):
yy = [1-u, 0, 1]
yy.sort()
m = yy[1]
# m = clamp(1-u, 0, 1);
m2 = m*m
return m2*m2*m # pow(m,5)
def sqrt(x):
return math.sqrt(x)
def log(x):
return math.log(x)
H = 29
S = np.random.uniform(0,1,1000000)
V = np.random.uniform(0,1,1000000)
metalic = np.random.uniform(0,1,1000000)
subsurface = np.random.uniform(0,1,1000000)
specular = np.random.uniform(0,1,1000000)
specularTint = np.random.uniform(0,1,1000000)
roughness = np.random.uniform(0,1,1000000)
anistrophic = np.random.uniform(0,1,1000000)
sheen = np.random.uniform(0,1,1000000)
sheenTint = np.random.uniform(0,1,1000000)
clearCoat = np.random.uniform(0,1,1000000)
clearCoatGloss = np.random.uniform(0,1,1000000)
result_all = []
ran = 0.00140
for i in range(1000000):
HSV2RGB(H,S[i],V[i])
baseColor0 = R
baseColor1 = G
baseColor2 = B
ax_ = max(.001, (roughness[i])**2/math.sqrt(1-anistrophic[i]*.9))
ay_ = max(.001, (roughness[i])**2*math.sqrt(1-anistrophic[i]*.9))
Cdlum_ = 0.3*baseColor0**2.2 + 0.6*baseColor1**2.2 + 0.1*baseColor2**2.2
BSDF0 = f0(baseColor0,baseColor1,baseColor2,metalic[i],subsurface[i],specular[i],specularTint[i], roughness[i],anistrophic[i],sheen[i],sheenTint[i],clearCoat[i],clearCoatGloss[i],ax_,ay_,Cdlum_)
BSDF1 = f1(baseColor0,baseColor1,baseColor2,metalic[i],subsurface[i],specular[i],specularTint[i], roughness[i],anistrophic[i],sheen[i],sheenTint[i],clearCoat[i],clearCoatGloss[i],ax_,ay_,Cdlum_)
BSDF2 = f2(baseColor0,baseColor1,baseColor2,metalic[i],subsurface[i],specular[i],specularTint[i], roughness[i],anistrophic[i],sheen[i],sheenTint[i],clearCoat[i],clearCoatGloss[i],ax_,ay_,Cdlum_)
BSDF3 = f3(baseColor0,baseColor1,baseColor2,metalic[i],subsurface[i],specular[i],specularTint[i], roughness[i],anistrophic[i],sheen[i],sheenTint[i],clearCoat[i],clearCoatGloss[i],ax_,ay_,Cdlum_)
BSDF4 = f4(baseColor0,baseColor1,baseColor2,metalic[i],subsurface[i],specular[i],specularTint[i], roughness[i],anistrophic[i],sheen[i],sheenTint[i],clearCoat[i],clearCoatGloss[i],ax_,ay_,Cdlum_)
BSDF5 = f5(baseColor0,baseColor1,baseColor2,metalic[i],subsurface[i],specular[i],specularTint[i], roughness[i],anistrophic[i],sheen[i],sheenTint[i],clearCoat[i],clearCoatGloss[i],ax_,ay_,Cdlum_)
BSDF6 = f6(baseColor0,baseColor1,baseColor2,metalic[i],subsurface[i],specular[i],specularTint[i], roughness[i],anistrophic[i],sheen[i],sheenTint[i],clearCoat[i],clearCoatGloss[i],ax_,ay_,Cdlum_)
BSDF7 = f7(baseColor0,baseColor1,baseColor2,metalic[i],subsurface[i],specular[i],specularTint[i], roughness[i],anistrophic[i],sheen[i],sheenTint[i],clearCoat[i],clearCoatGloss[i],ax_,ay_,Cdlum_)
BSDF8 = f8(baseColor0,baseColor1,baseColor2,metalic[i],subsurface[i],specular[i],specularTint[i], roughness[i],anistrophic[i],sheen[i],sheenTint[i],clearCoat[i],clearCoatGloss[i],ax_,ay_,Cdlum_)
BSDF9 = f9(baseColor0,baseColor1,baseColor2,metalic[i],subsurface[i],specular[i],specularTint[i], roughness[i],anistrophic[i],sheen[i],sheenTint[i],clearCoat[i],clearCoatGloss[i],ax_,ay_,Cdlum_)
BSDF10 = f10(baseColor0,baseColor1,baseColor2,metalic[i],subsurface[i],specular[i],specularTint[i], roughness[i],anistrophic[i],sheen[i],sheenTint[i],clearCoat[i],clearCoatGloss[i],ax_,ay_,Cdlum_)
BSDF11 = f11(baseColor0,baseColor1,baseColor2,metalic[i],subsurface[i],specular[i],specularTint[i], roughness[i],anistrophic[i],sheen[i],sheenTint[i],clearCoat[i],clearCoatGloss[i],ax_,ay_,Cdlum_)
if (BSDF0 > -ran) and (BSDF0 < ran) and (BSDF1 < ran) and (BSDF1 > -ran) and (BSDF2 < ran) and (BSDF2 > -ran) and (BSDF3 < ran) and (BSDF3 > -ran) and (BSDF4 < ran) and (BSDF4 < ran) and (BSDF5 > -ran) and (BSDF5 < ran) and (BSDF6 > -ran) and (BSDF6 < ran) and (BSDF7 > -ran) and (BSDF7 < ran) and (BSDF8 < ran) and (BSDF8 > -ran) and (BSDF9 < ran) and (BSDF10 < ran) and (BSDF10 > -ran) and (BSDF11 < ran):
result = []
result.append(baseColor0)
result.append(baseColor1)
result.append(baseColor2)
result.append(metalic[i])
result.append(subsurface[i])
result.append(specular[i])
result.append(specularTint[i])
result.append(roughness[i])
result.append(anistrophic[i])
result.append(sheen[i])
result.append(sheenTint[i])
result.append(clearCoat[i])
result.append(clearCoatGloss[i])
result_all.append(result)
else:
pass
print(len(result_all))
print(result_all)