通过MERL100计算Blender Disney BRDF参数

本文详细介绍了如何使用MERL100数据库来计算Blender中DisneyPrincipledBSDF的材质参数,以实现更真实的渲染效果。作者通过分析数据库、编译C++程序、调整色彩校正参数,以及估算光源系数,最终得出适用于特定材质的参数。整个过程涉及随机数解方程、优化算法等技术,以求得与实际观察结果相匹配的材质表现。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

图:大理石参数渲染结果

概述

我在刚开始接触十多个参数的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-235cherry235(木材清漆木材
gold-metallic-paint2金色金属漆2涂料colonial-maple-223CM223(木材清漆木材
gold-metallic-paint3金色金属漆3涂料fruitwood-241fw241(木材清漆木材
gold-paint金色涂料涂料ipswich-pine-221ip221(木材清漆木材
green-latex绿色乳胶涂料natural-209n209(木材清漆木材
green-metallic-paint绿色金属涂料涂料pickled-oak-260po260(木材清漆木材
green-metallic-paint2绿色金属涂料2涂料special-walnut-224sw224(木材清漆木材
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灰色塑料塑料ss440ss440碳素钢金属
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)

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值