CN104700445B

CN104700445B - A kind of BRDF reflection model deriving methods based on measurement data

Info

Publication number: CN104700445B
Application number: CN201510128193.3A
Authority: CN
Inventors: 王璐; 李冰辉; 徐延宁; 孟祥旭; 陆巧; 孟祥雨
Original assignee: Shandong University
Current assignee: Shandong University
Priority date: 2015-03-23
Filing date: 2015-03-23
Publication date: 2018-07-27
Anticipated expiration: 2035-03-23
Also published as: CN104700445A

Abstract

The invention discloses a kind of BRDF reflection model deriving methods based on measurement data, include the following steps：Step 1, Data Dimensionality Reduction：Based on BRDF databases, high dimensional data is subjected to dimensionality reduction, indicates that high dimension vector, BRDF are writing a Chinese character in simplified form for bidirectional reflectance function with low-dimensional vector；Step 2, the calculating of eigen vector：The object properties estimated performance direction that user is specified with support vector machine method；Step 3, derivative new BRDF：The characteristic direction calculated in step 2 is applied on the initial BRDF selected in step 1, new BRDF is derived；Step 4 is rendered using derivative new BRDF.The result of rendering is very true to nature.The present invention allows a user to specify that a series of intuitive parameters are used for changing the attribute of BRDF.The parameter with perception meaning used in the present invention is easier to control and uses.The present invention need not explicitly store all BRDF measurement data, save memory space.

Description

A kind of BRDF reflection model deriving methods based on measurement data

Technical field

The present invention relates to graphics Realistic Rendering fields, and in particular to a kind of bidirectional reflectance function based on measurement data (BRDF) reflection model deriving method includes dimensionality reduction, reconstruction, derivative and rendering side based on original BRDF measurement data Method.

Background technology

In general, being provided using the analysis reflection model (such as Cook-Torrance) that physics inspires in computer graphics BRDF, these BRDF models are the approximation of the reflection of real-world object.Most of analysis reflection models are only limitted to describe specific The material of type, a given model are only capable of indicating the phenomenon designed by it.The parameter of this model can pass through in principle True measure obtains, but is but difficult to carry out in practice.

Another measurement model method getparms is directly to be obtained from true BRDF samplings, is then used various Optimisation technique carries out data fitting using certain analysis reflection model.This method being fitted afterwards that first measures has several drawbacks in that.It is first First, the approximation only really reflected by the BRDF calculated with analytic function, the BRDF values of measurement are simultaneously imprecise equal to analysis The value of model.It measures again approximating method and assumes that there are noises, the process of fitting to have filtered out these noises during measurement. For in terms of this, this will neglect many modeling errors caused by approximate analysis reflection model.Many objects it is anti- Penetrating characteristic may be within the scope of these modeling error.Secondly, the performance that Select Error function optimizes on the mold is simultaneously It is undesirable.Finally, it cannot be guaranteed that the result of optimization can generate best model.Since most of BRDF models are nonlinearities , the Optimization Framework used in fit procedure is largely dependent upon the initial guess of model parameter, these are initially guessed The quality of survey will have significant impact to model final argument value.

Another method is to obtain intensive BRDF to measure, and is directly rendered using these measurement data.This side Method remains the measurement data lost in data fitting method.However, this process takes very much, because it is needed to scene In all objects carry out reflection measurement.In addition, the change of any thingness, will need to pick up one and expectation category The consistent object of property, then obtains its reflection measurement.

In existing patent, patent disclosed in 2012, it is proposed that a kind of indoor full automatic BRDF measuring devices, including branch Tablet is supportted, has horizontal revolving stage, sample stage to be mounted in the central through hole of horizontal revolving stage in support platform, mirror is installed on horizontal revolving stage Head bracket installs camera lens pivoted arm on lens bracket, and spectrometer camera lens is installed at camera lens pivoted arm arm end, and the optical axis of spectrometer camera lens is always It is directed toward the sample stage center, light source guide rail is also installed in support platform, artificial light source, artificial light source light are set on light source guide rail Axis is directed toward sample stage center always.The BRDF measuring devices can effectively improve positioning accuracy, reduce measurement period, same to time Source collimation and uniformity are improved, and then are effectively guaranteed measurement accuracy.Since the present apparatus can only acquire in reality Some objects, so using above there is significant limitation.

In conclusion specific defect of the existing technology is：

1. the material parameter that existing model uses, however, it would be possible to measure, but be difficult to obtain in practice.

2. measuring approximating method assumes that there are noises, the process of fitting to have filtered out these noises in measurement process, this Many modeling errors as caused by approximate analysis reflection model will be neglected.

3. obtaining the method that intensive BRDF is measured, this process takes very much, in addition, the applicability of this method is not By force, the change of any thingness will need to pick up an object consistent with desired properties, then obtain the anti-of it Penetrate rate.

Invention content

To solve the shortcomings of the prior art, the invention discloses a kind of application method based on BRDF measurement data, Including BRDF Data Dimensionality Reductions, reconstruction and deriving method, it can be used for the rendering of high presence material.Whole process of the present invention In two steps：The measurement of BRDF and the use of BRDF.Here the use of BRDF is mainly introduced, part is measured and is mainly managed using Massachusetts Measurement data disclosed in engineering college.

To achieve the above object, concrete scheme of the invention is as follows：

A kind of BRDF reflection model deriving methods based on measurement data, include the following steps：

Step 1, Data Dimensionality Reduction：Based on BRDF databases, high dimensional data is subjected to dimensionality reduction, higher-dimension is indicated with low-dimensional vector Vector, BRDF are writing a Chinese character in simplified form for bidirectional reflectance function；

Step 2, the calculating of eigen vector：The object properties estimated performance side that user is specified with support vector machine method To；

Step 3, derivative new BRDF：By the characteristic direction calculated in step 2 apply in step 1 select it is initial On BRDF, new BRDF is derived；

Step 4 is rendered using derivative new BRDF.

In the step 1, the mode that high dimensional data carries out dimensionality reduction uses linear mode or nonlinear way, using linear When mode dimensionality reduction, PCA, that is, principal component analytical method is used.

When the use linear mode dimensionality reduction, using PCA, that is, principal component analytical method, specially：

U=[u1, u2 ..., uD] is enabled to indicate new coordinate system, then X=[x1, x2 .., xN], N indicate BRDF files Number, that is, the number put, then enabling Y=U*X, because U is symmetrical matrix, U'*U=I, matrix X can pass through X= U'*Y is calculated；Due to indicating that D is tieed up with d dimensions, so only indicating former vector with the wherein d of base vector；

Wherein, y is enabled_iIndicate a point under D dimension spaces, wherein u_jJ-th of base vector under denotation coordination system U, x_jiIt indicates The j-th point of projection under coordinate system on i-th of base vector direction, d indicate the number of new coordinate system base vector,It indicates only The y indicated is only mapped with d base vector_i, Y indicate it is N number of point composition matrix.

When estimation is using linear mode dimensionality reduction, the information loss of the entire data set of PCA, that is, principal component analytical method is used：

Formula (2)

Wherein, ∑ is the covariance matrix of all the points in data set, and N indicates the number of point, E_MIt indicates information loss, uses T To indicate the transposed matrix of some matrix.

Work as u_jWhen the characteristic value obtained corresponding to covariance feature matrix ∑ Eigenvalues Decomposition, data degradation reaches minimum：

U Λ=∑ U formula (5)

Λ=U^T∑ U formula (6)

Wherein, Λ is the eigenvalue λ successively decreased_jDiagonal matrix, U is the orthogonal moment that the corresponding feature vector of characteristic value indicates Battle array；The information loss of data set can be reduced to following form：

When d feature vector being chosen is the maximum characteristic values of preceding d corresponding feature vector, data degradation reaches most It is small.When D reaches million grades of sizes, the Eigenvalues Decomposition of matrix B and it is defined first：

B=Y^TY=V Λ_BV^TFormula (8)

Wherein, V is the orthogonal matrix of feature vector composition, Λ_BThe eigenvalue cluster that successively decreases at diagonal matrix, then, One Y of premultiplication, the right side multiply a V, then

YY_TYV=YV Λ_BV^TV formula (9)

∑ (YV)=(YV) Λ_BFormula (10)

It converts again,

And because

So

Λ=Λ_BFormula (15).

In the step 2, the calculating of eigen vector is entirely counting not in current BRDF come according to some characteristic Sorted out according to collection, if if, it being labeled as 1, -1 is labeled as not if；Algorithm of support vector machine is found two Data set gives separated optimal hyperlane.

In step 2, the calculating of eigen vector, specially：

First, it whether there is in some BRDF according to specified characteristic, entire data set can be classified as two classes：

(x₁,y₁),...,(x_n,y_n),x_i∈R^d,y∈{-1,1}

Assuming that the two classification are linear separabilities, then a hyperplane may be found by data set to separately：

y_i[w^Tx_i+ b] >=1, i=1 ..., n.

Wherein w^Tx_i+ b means that a hyperplane, and the distance of nearest point is equal in hyperplane to two datasets；

Wherein, w indicates that the normal vector of hyperplane, b represent a constant, x_iIndicate at i-th point, w^TIndicate the transposition of w.Meter When calculating w and b：Meet：

subject to y_i(w^Tx_i+ b) >=1, i=1 ..., n.

The optimal value of w and b, which calculates, to be calculated by introducing Lagrange multiplier：

subject to α_i>=0, i=1 ..., n.

Wherein, α indicates the vector that n Lagrange multiplier is constituted, α_iIndicate that i-th of Lagrange multiplier, J indicate to use Parameter w, the object function that b, α are indicated.

Local derviation is asked to object function：

Then, the problem of optimization, can be deformed into following form：

α_i>=0, i=1 ..., n,

The optimization of this problem can be solved by quadratic programming.

Quadratic programming is solved, and specifically used SMO algorithms, problem are again converted to following form, KKT conditions：

The case where for (a), shows α_iIt is normally to classify, in border inner；

The case where for (b), shows α_iIt is supporting vector, on boundary；

The case where for (c), shows α_iIt is between two boundaries；

Optimal solution needs to meet KKT conditions, i.e. (a) (b) (c) three conditions will meet.

Find out the α for being unsatisfactory for KKT conditions_iAnd be updated, but α_iAlso another is constrained, i.e.,

Wherein l indicates the number of point.

The present invention updates α by another method_iAnd α_j：

Ensure that in this way and for 0 constraint.Whereinα before indicating to update_iValue,Before indicating to update α_jValue,Indicate α after updating_iValue,Indicate α after updating_jValue.

Utilize y_ia_i+y_ja_j=constant eliminates α_i, one can be obtained about single argument α_jA convex quadratic programming ask Topic does not consider that it constrains 0≤α_j≤ C, the solution for obtaining it are：

Wherein, C indicates penalty coefficient, i.e., if some point belongs to certain one kind, but it deviates from such, goes on boundary Or the place of other classes is gone, C shows more to be not desired to abandon this point more greatly, and boundary will reduce.It indicates after updating α_jValue, α_jIndicate the value before update, y_iIndicate belong to which kind of at i-th point.

In addition,

E_i=u_i-y_i

η=k (x_i,x_i)+k(x_j,x_j)-2k(x_i,x_j)

Consider 0≤α_jAfter the constraint of≤C, the solution for obtaining it is

Wherein,Expression applies α after constraint_jValue,

So,

Find such a pair of α_iAnd α_j：

α_iThe multiplier of KKT conditions, α can be unsatisfactory for by looking for_jCan look for and meet condition max | E_i-E_j| multiplier.

The update of b：

Wherein,

According to the method, all Lagrange multiplier and b have been updated, you can calculate final eigen vector：

W=∑s α_iy_ix_i。

In the step 4, is rendered, specifically included using derivative BRDF：

Step 4.1 directly inquires the brdf of BRDF file acquisition corresponding points, shown herein as the reflection under current light source angle Rate information calculates the bloom of corresponding points, that is, diffuse=Cl*brdf*L.N further according to Lighting information.Wherein, Cl generations The color in mass color source, L represent the direction vector of light source, and N represents the current normal direction for rendering point.

Step 4.2, the reflection that object is calculated by BRDF is calculated using importance sampling and is reflected.

The illumination of current point is other than direct illumination, and also other shine its indirect light, so, product can be passed through Point method calculate the illumination of current point：

L_r(θ_r,φ_r)=∫_Ωf_r(θ_h,θ_d,φ_d)L_i(θ_i,φ_i)cosθ_idω_iFormula (16)

Wherein, Ω defines incident light hemisphere, L_rIt is the irradiation level on expected reflection direction, ω represents all directions, θ_hTable Show angle, the θ of h and normal direction_dRepresent θ_iWith θ_rDifference, θ_iAnd θ_rRespectively direction of visual lines and normal direction, incident light direction and normal direction Angle, φ_dIndicate φ_iWith φ_rDifference, φ_iAnd φ_rThe respectively projection of direction of visual lines in the horizontal plane and t axis and incidence The angle of the projection and t axis of light direction in the horizontal plane, f_rIt is reflective function, with brdf described above, it can be according to ginseng Number θ_h, θ_d, φ_dIt tables look-up to obtain.

According to the thought of Monte Carlo, principal value of integral can be estimated by sampling N number of point, the value of N is bigger, the knot of estimation Fruit is more accurate：

Wherein p (x) represents probability density function, and f (x) indicates the integrand of Monte Carlo integral, F_NIndicate estimation Integrated value.Density function is unknown, when calculating density function：

Formula (18) is the optimal selection of density function, and f (x) is integrand, so formula may be used in the present invention (19) density function as the present invention, wherein in formula (19), θ_iAnd θ_rRespectively direction of visual lines and normal direction, incident light direction With the angle of normal direction, φ_diffIndicate the difference of the φ values of exit direction and the φ values of incident direction, i.e. φ_r-φ_i, f_sum(θ_r,θ_i, φ_diff) be tri- components of corresponding points BRDF and.

θ_iAnd φ_diffAll it is unknown, the first step, calculates θ_i, second step, the θ found out according to back_iCalculate φ_diff：

Calculate θ_i：

In this step, present invention uses the Marginal density function,s of formula (21)：

In order to keep sampling process more efficient, the present invention needs to store per a pair of θ_rAnd θ_iCorresponding edge density value；

The function is a Marginal density function, in order to choose θ according to importance_i, the PDF for further calculating out it is Probability density function, each value of probability density function are the numerical value between 0 to 1, for the ease of coming according to its size Decide whether to choose corresponding θ_i, CDF i.e. cumulative density function is further calculated, then, so that it may to generate between one 0 to 1 Random number, probability is big, and the probability arrived at random is also big.

Calculate φ_diff：

In this step, P (φ are calculated with formula (22)_diff|θ_r,θ_i), CDF is further calculated, in order to efficiently consider, is needed Storage is per a pair of θ_r, θ_iCorresponding CDF generates the random number between one 0 to 1, then uses the method for binary chop in table Find out corresponding φ_diff；

So far, θ_iAnd φ_diffIt all calculates, substitutes into formula (19), calculate corresponding probability；Finally substitute into Inside formula (18), that is, reflection results are acquired, by the direct illumination results added of reflection results and step 4.1, you can wrapped Global illumination result containing reflective information.

Beneficial effects of the present invention：

The present invention uses the BRDF measured databases that MIT is provided, by Data Dimensionality Reduction, reconstruction, finally deriving new BRDF finally applies to the BRDF in Realistic Rendering scene.The invention has the advantage that：

1. the result rendered is very true to nature.

2. the present invention allows a user to specify that a series of intuitive parameters are used for changing the attribute of BRDF.

3. the parameter with perception meaning used in the present invention is easier to control and uses.

4. the present invention need not explicitly store all BRDF measurement data, memory space is saved.

Description of the drawings

Fig. 1 is the definition figure of BRDF；

Fig. 2 a are conventional coordinates；

The coordinate system that Fig. 2 b change；

Fig. 3 be from left to right BRDF rendering effects using 5,10,20,30,45,60,100 main feature reconstructions respectively Figure.

Fig. 4 is the rendering effect figure that bloom effect is gradually added in the BRDF of the coarse effect of blue.

Wherein, in Fig. 1, L_iFor the direction of realization, L_rFor the direction of incident light, θ_iAnd θ_rRespectively direction of visual lines and Z axis The angle of angle and incident light direction and Z axis, φ_iAnd φ_rThe respectively angle of direction of visual lines and X-axis and incident light direction and X-axis Angle.

In Fig. 2 a, ω_iAnd ω_rRespectively direction of visual lines and incident light direction, θ_iAnd θ_rRespectively direction of visual lines and normal direction, The angle of incident light direction and normal direction, φ_iAnd φ_rThe respectively projection of direction of visual lines in the horizontal plane and t axis and incident light direction The angle of projection and t axis in the horizontal plane, f_rIt is reflective function, it can be according to parameter θ_h, θ_d, φ_dIt tables look-up to obtain.

In Fig. 2 b, ω_iAnd ω_rMeaning is identical as left figure, and h represents ω_iAnd ω_rAngular bisector, θ_hAnd φ_hRepresent h and method To the projection in the horizontal plane of angle and h and the angle of t axis, θ_dAnd φ_dRepresent θ_iWith θ_rDifference and φ_iWith φ_rDifference.

Specific implementation mode：

The present invention is described in detail below in conjunction with the accompanying drawings：

Step 1. Data Dimensionality Reduction.

Linear mode and nonlinear way can be used to the dimensionality reduction of high dimensional data, herein, the present invention discusses linear side Formula dimensionality reduction.Linear mode dimensionality reduction at present using it is more be PCA (principal component analysis) method.

The main target of PCA is to find a linear transformation, data from High Dimensional Mapping to low-dimensional, while preserving number again According to most of information.Yi is enabled to indicate a point under D dimension spaces, then 100 BRDF files can form a D*100 Matrix.U=[u1, u2 ..., uD] is enabled to indicate new coordinate system, then X=[x1, x2 .., xn], n indicate BRDF files Number, that is, the number put.So Y=U*X.Because U is symmetrical matrix, U'*U=I.Matrix X can pass through X= U'*Y is calculated.

Due to indicating that D is tieed up with d dimensions, so only indicating former vector with the wherein d of base vector.

The information loss of entire data set can be estimated with the following formula：

Wherein, ∑ is the covariance matrix of all the points in data set, and N indicates the number of point, E_MIndicate information loss.Use T To indicate the transposed matrix of some matrix.

Work as u_jWhen the characteristic value obtained corresponding to covariance feature matrix ∑ Eigenvalues Decomposition, data degradation reaches minimum.

U Λ=∑ U formula (5)

Λ=U^T∑ U formula (6)

Wherein, Λ is the eigenvalue λ successively decreased_jDiagonal matrix, U is the orthogonal moment that the corresponding feature vector of characteristic value indicates Battle array.

The information loss of data set can be reduced to following form：

This also means that, d feature vector being chosen of the present invention be the corresponding feature of preceding d maximum characteristic values to Amount can thus make information loss reach minimum.

However, when D is prodigious (a BRDF file in this program can reach 180*90*90*3 dimensions), it is very The difficult Eigenvalues Decomposition that can not possibly calculate covariance matrix and it in other words.

So the present invention alternatively calculates indirectly.The present invention defines matrix B and its characteristic value point first Solution.

B=Y^TY=V Λ_BV^TFormula (8)

Wherein, V is the orthogonal matrix of feature vector composition, Λ_BThe eigenvalue cluster that successively decreases at diagonal matrix.Then, One Y of premultiplication, the right side multiply a V, then

YY_TYV=YV Λ_BV^TV formula (9)

∑ (YV)=(YV) Λ_BFormula (10)

It converts again,

And because

So

Λ=Λ_BFormula (15)

When the present invention realizes, characteristic value and spy are calculated using SVD (singular value decomposition) methods of matlab Sign vector.

The calculating of step 2. eigen vector

The present invention needs sorting out entire data set not in current BRDF according to some characteristic, if If, it is labeled as 1, -1 is labeled as not if, then can be carried out the calculating of eigen vector.

Algorithm of support vector machine can be found two datasets to separated optimal hyperlane.

First, it whether there is in some BRDF according to specified characteristic, entire data set can be classified as two classes.

(x₁,y₁),...,(x_n,y_n),x_i∈R^d,y∈{-1,1}

By taking green attribute as an example, 100 BRDF files can be classified as two classes.

[-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,- 1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,1,1,1,1,1,-1,-1,-1,-1,-1,-1,-1,-1,-1,- 1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1]；

Wherein 1 to represent the BRDF include green attribute, and -1, which represents the BRDF, does not include green attribute, or perhaps green Attribute unobvious.

Assuming that the two classification are linear separabilities, then a hyperplane may be found by data set to separately.

y_i[w^Tx_i+ b] >=1, i=1 ..., n.

Wherein w^Tx_i+ b means that a hyperplane.The distance of nearest point is equal in hyperplane to two datasets.

In order to calculate w and b, need to solve the problems, such as following：

subject to y_i(w^Tx_i+ b) >=1, i=1 ..., n.

The optimal value of w and b, which calculates, to be calculated by introducing Lagrange multiplier.

subject to α_i>=0, i=1 ..., n.

Local derviation is asked to object function：

Then, the problem of optimization, can be deformed into following form：

α_i>=0, i=1 ..., n,

The optimization of this problem can be solved by quadratic programming.Herein, present invention uses SMO algorithms.

Problem is again converted to following form (KKT conditions)：

The case where for (a), shows α_iIt is normally to classify, in border inner (the point y correctly to classify_i*f(i)≥0)

The case where for (b), shows α_iIt is supporting vector, on boundary.

The case where for (c), shows α_iIt is between two boundaries.

So to find out the α for being unsatisfactory for KKT conditions_iAnd it is updated.But α_iAlso another is constrained, i.e.,

Wherein l indicates the number of point.

The present invention updates α by another method_iAnd α_j：

Ensure that in this way and for 0 constraint.Whereinα before indicating to update_iValue,α before indicating to update_j Value,Indicate α after updating_iValue,Indicate α after updating_jValue.

Wherein, C indicates penalty coefficient, i.e., if some point belongs to certain one kind, but it deviates from such, goes on boundary Or the place of other classes is gone, C shows more to be not desired to abandon this point more greatly, and boundary will reduce.It indicates after updating α_jValue, α_jIndicate the value before update, y_iIndicate belong to which kind of at i-th point.In addition,

E_i=u_i-y_i

η=k (x_i,x_i)+k(x_j,x_j)-2k(x_i,x_j)

Consider 0≤α_jAfter the constraint of≤C, the solution for obtaining it is

Wherein,Expression applies α after constraint_jValue,

So,

Find such a pair of α_iAnd α_j：

The update of b：

Wherein,

W=∑s α_iy_ix_i

W=[- 0.0003,0.0006,0.0003,0.0007,0.0000,0.0004,0.0010,0.0000 ,- 0.0007,-0.0002,0.0000,0.0003,-0.0002,0.0001,0.0002,-0.0006,-0.0001,0.0002, 0.0002,-0.0001,0.0001,0.0004,0.0002,-0.0002,0.0000,-0.0001,-0.0001,0.0000- 0.0001,0.0000,-0.0000,0.0001,0.0000,0.0000,-0.0003,-0.0000,-0.0002-0.0000,- 0.0001,0.0000,0.0000,-0.0000,-0.0000,0.0001,-0.0000]

The eigen vector is multiplied by certain coefficient by step 3., in the BRDF space vectors after the dimensionality reduction that is added to, then is rebuild To luv space, the new BRDF that can directly use just has been obtained, will show green category with new BRDF files rendering Property.

Step 4. is rendered using derivative BRDF.

Step 4.2, the reflection of object is calculated by BRDF, hereafter main introduce calculates reflection using importance sampling.

L_r(θ_r,φ_r)=∫_Ωf_r(θ_h,θ_d,φ_d)L_i(θ_i,φ_i)cosθ_idω_iFormula (16)

The illumination of current point is other than direct illumination, and also other shine its indirect light, so, product can be passed through Point method calculate the illumination of current point, such as formula (16), wherein Ω defines incident light hemisphere, L_rIt is in expected reflection direction On irradiation level, ω represents all directions, θ_h、θ_d、φ_dIntroduction please see Figure 2 introduction.

According to the thought of Monte Carlo, principal value of integral can be estimated by sampling N number of point, the value of N is bigger, the knot of estimation Fruit is more accurate, and such as formula (17), wherein p (x) represents probability density function.

Density function is unknown, introduces how to calculate density function below.

Formula (18) is the optimal selection of density function, and f (x) is integrand, so formula may be used in the present invention (19) density function as the present invention, the wherein θ in formula (19)_rAnd θ_iIt is identical as the meaning in Fig. 2, φ_diffWith in Fig. 2 φ_dMeaning is identical.

f_sum(θ_r,θ_i,φ_diff) be tri- components of corresponding points BRDF and.

Now, θ_iAnd φ_diffAll it is unknown, so the present invention can calculate in two steps.The first step calculates θ_i, second Step, the θ found out according to back_iCalculate φ_diff。

Step 4.2.1 calculates θ_i

In order to keep sampling process more efficient, the present invention needs to store per a pair of θ_rAnd θ_iCorresponding edge density value.

Note that the function is a Marginal density function, in order to choose θ according to importance_i, the present invention also needs into one Step calculates its PDF (probability density function).Each value of probability density function is the numerical value between 0 to 1, in order to just In deciding whether to choose corresponding θ according to its size_i, the present invention needs to further calculate CDF (cumulative density function).So Afterwards, the present invention can generate the random number between one 0 to 1, and probability is big, and the probability arrived at random is also big, and here it is importance Where the marrow of sampling.

Step 4.2.2 calculates φ_diff

In this step, P (φ are calculated with formula (22)_diff|θ_r,θ_i), CDF is further calculated, in order to efficiently consider, is needed Storage is per a pair of θ_r, θ_iCorresponding CDF generates the random number between one 0 to 1, then uses the method for binary chop in table Find out corresponding φ_diff。

The result that the present invention renders is more life-like compared to traditional parsing reflection model.Method used in the present invention All BRDF measurement data need not be stored, compared to for traditional measurement storage method, save memory space.The present invention State modulator it is more intuitive, easier, while a large amount of BRDF can also be derived for users to use.

Above-mentioned, although the foregoing specific embodiments of the present invention is described with reference to the accompanying drawings, not protects model to the present invention The limitation enclosed, those skilled in the art should understand that, based on the technical solutions of the present invention, those skilled in the art are not Need to make the creative labor the various modifications or changes that can be made still within protection scope of the present invention.

Claims

1. a kind of BRDF reflection model deriving methods based on measurement data, characterized in that include the following steps：

Step 1, Data Dimensionality Reduction：Based on BRDF databases, high dimensional data is subjected to dimensionality reduction, high dimension vector is indicated with low-dimensional vector, BRDF is writing a Chinese character in simplified form for bidirectional reflectance function；

Step 2, the calculating of eigen vector：The object properties estimated performance direction that user is specified with support vector machine method；

Step 3, derivative new BRDF：The characteristic direction calculated in step 2 is applied to the initial BRDF selected in step 1 On, derive new BRDF；

Step 4 is rendered using derivative new BRDF；

In the step 4, is rendered, specifically included using derivative BRDF：

Step 4.1 directly inquires the brdf of BRDF file acquisition corresponding points, believes shown herein as the reflectivity under current light source angle Breath calculates the bloom of corresponding points, that is, diffuse=Cl*brdf*L.N, wherein Cl represents light further according to Lighting information The color in source, L represent the direction vector of light source, and N represents the current normal direction for rendering point；

Step 4.2, the reflection that object is calculated by BRDF is calculated using importance sampling and is reflected, by reflection results and step 4.1 Direct illumination results added, you can obtain including the global illumination result of reflective information；

The illumination of current point is other than direct illumination, and also other shine its indirect light, so, pass through the method for integral To calculate the illumination of current point：

L_r(θ_r,φ_r)=∫_Ωf_r(θ_h,θ_d,φ_d)L_i(θ_i,φ_i)cosθ_idω_iFormula (16)

Wherein, Ω defines incident light hemisphere, L_rIt is the irradiation level on expected reflection direction, ω represents all directions, θ_hRepresent h With the angle of normal direction, θ_dAnd φ_dRepresent θ_iWith θ_rDifference and φ_iWith φ_rDifference, f_rIt is reflective function；θ_iWith θ_rRespectively regard Line direction and normal direction, the angle of incident light direction and normal direction；φ_iWith φ_rRespectively direction of visual lines projection in the horizontal plane and t The angle of the projection and t axis of axis and incident light direction in the horizontal plane.

2. a kind of BRDF reflection model deriving methods based on measurement data as described in claim 1, characterized in that the step In rapid one, the mode that high dimensional data carries out dimensionality reduction uses linear mode or nonlinear way, when using linear mode dimensionality reduction, uses PCA, that is, principal component analytical method.

3. a kind of BRDF reflection model deriving methods based on measurement data as claimed in claim 2, characterized in that described to adopt When with linear mode dimensionality reduction, using PCA, that is, principal component analytical method, specially：U=[u1, u2 ..., uD] is enabled to indicate new seat Mark system, then X=[x1, x2 .., xn], n indicate the number of BRDF files, that is, the number put, then Y=U*X, because of U It is symmetrical matrix, so U'*U=I, matrix X is calculated by X=U'*Y；Due to indicating that D is tieed up with d dimensions, so only using basal orientation The wherein d of amount indicates former vector；

Wherein, y is enabled_iIndicate a point under D dimension spaces, wherein u_jJ-th of base vector under denotation coordination system U, x_jiIt indicates j-th Projection of the point under coordinate system on i-th of base vector direction, d indicate the number of new coordinate system base vector,D is only used in expression The y that a base vector mapping indicates_i, Y indicate it is N number of point composition matrix.

4. a kind of BRDF reflection model deriving methods based on measurement data as described in claim 1, characterized in that the step In rapid two, the calculating of eigen vector is being sorted out entire data set according to some characteristic not in current BRDF, such as Fruit is labeled as 1 if, it, and -1 is labeled as not if；Algorithm of support vector machine is found two datasets to separated Optimal hyperlane.

5. a kind of BRDF reflection model deriving methods based on measurement data as claimed in claim 4, characterized in that step 2 In, the calculating of eigen vector, specially：

First, it will be whether there is in some BRDF according to specified characteristic, entire data set be classified as two classes：

(x₁,y₁),...,(x_n,y_n),x_i∈R^d,y∈{-1,1}

y_i[w^Tx_i+ b] >=1, i=1 ..., n.

Wherein, w indicates that the normal vector of hyperplane, b represent a constant, x_iIndicate at i-th point, w^TIndicate the transposition of w.

6. a kind of BRDF reflection model deriving methods based on measurement data as described in claim 1, characterized in that according to illiteracy The thought of special Carlow estimates principal value of integral by sampling N number of point, and the value of N is bigger, and the result of estimation is more accurate：

Wherein p (x) represents probability density function, and f (x) indicates the integrand of Monte Carlo integral, F_NIndicate the integral of estimation Value.

7. a kind of BRDF reflection model deriving methods based on measurement data as claimed in claim 6, characterized in that density letter Number is unknown, when calculating density function：

Formula (18) is the optimal selection of density function, and f (x) is integrand, and density function is used as using formula (19), wherein θ in formula (19)_iAnd θ_rRespectively direction of visual lines and normal direction, the angle of incident light direction and normal direction, φ_diffIndicate outgoing side To φ values and incident direction φ values difference, i.e. φ_r-φ_i, f_sum(θ_r,θ_i,φ_diff) it is tri- components of corresponding points BRDF With.

8. a kind of BRDF reflection model deriving methods based on measurement data as claimed in claim 7, characterized in that θ_iWith φ_diffAll it is unknown, the first step, calculates θ_i, second step, the θ found out according to back_iCalculate φ_diff：

Calculate θ_i：

Storage is per a pair of θ_rAnd θ_iCorresponding edge density value；The Marginal density function, of formula (21) is used；

The function is a Marginal density function, in order to choose θ according to importance_i, further calculate out its PDF, that is, probability Density function, each value of probability density function are the numerical value between 0 to 1, for the ease of being determined according to its size Whether corresponding θ is chosen_i, further calculate CDF i.e. cumulative density function, then, so that it may with generate one 0 to 1 between with Machine number, probability is big, and the probability arrived at random is also big；

Calculate φ_diff：

P (φ are calculated with formula (22)_diff|θ_r,θ_i), CDF is further calculated, in order to efficiently consider, needs to store every a pair of θ_r, θ_i Corresponding CDF is generated the random number between one 0 to 1, is then found out in table using the method for binary chop corresponding φ_diff。