-
Notifications
You must be signed in to change notification settings - Fork 15
/
Copy pathdj_brdf.h
3673 lines (3203 loc) · 154 KB
/
dj_brdf.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/* dj_brdf.h - public domain BRDF toolkit
by Jonathan Dupuy
Do this:
#define DJ_BRDF_IMPLEMENTATION 1
before you include this file in *one* C++ file to create the implementation.
INTERFACING
define DJB_ASSERT(x) to avoid using assert.h.
define DJB_LOG(format, ...) to use your own logger (default prints in stdout)
QUICK NOTES
1. This is research code ;)
2. I have introduced sigma functions in the microfacet BRDF interface.
These functions are used for the Smith term because they are more
numerically stable than the G1 functions I defined in my thesis,
especially in noncentral configurations. This approach also allows
to avoid singulatities when either theta_o or theta_i approach pi/2.
In microflake theory, the sigma function corresponds to the projected
area of the flakes in a specific direction (see, e.g., "The SGGX
Microflake Distribution").
3. For people who are used to light transport / Mitsuba BSDF conventions:
I follow a "reversed" convention as:
=> the parameter 'i' denotes the direction towards the light source
=> the parameter 'o' denotes the direction towards the viewer
ACKNOWLEDGEMENTS
- Jiri Filip (provided code for djb::utia)
- Joel Kronander (provided parameters for djb::abc)
- Mahdi Bagher, Cyril Soler, Nicolas Holzschuch (provided parameters for djb::sgd)
*/
#ifndef DJB_INCLUDE_DJ_BRDF_H
#define DJB_INCLUDE_DJ_BRDF_H
#include <vector>
#include <string>
namespace djb {
/* Floating point precision */
#if DJB_USE_DOUBLE_PRECISION
typedef double float_t;
#else
typedef float float_t;
#endif
#ifndef DJB_EPSILON
#define DJB_EPSILON (float_t)1e-4
#endif
/* Exception API */
struct exc : public std::exception {
exc(const char *fmt, ...);
virtual ~exc() throw() {}
const char *what() const throw() {return m_str.c_str();}
std::string m_str;
};
/* Standalone vec3 utility */
struct vec3 {
static vec3 from_raw(const double *v) {return vec3((float_t)v[0], (float_t)v[1], (float_t)v[2]);}
static vec3 from_raw(const float *v) {return vec3((float_t)v[0], (float_t)v[1], (float_t)v[2]);}
static const float_t *to_raw(const vec3& v) {return &v.x;}
explicit vec3(float_t x = 0): x(x), y(x), z(x) {}
vec3(float_t x, float_t y, float_t z) : x(x), y(y), z(z) {}
explicit vec3(float_t theta, float_t phi);
float_t intensity() const {return (float_t)0.2126 * x + (float_t)0.7152 * y + (float_t)0.0722 * z;}
float_t x, y, z;
};
/* BRDF interface */
class brdf {
public:
// evaluate f_r
virtual vec3 eval(const vec3& i, const vec3& o,
const void *user_param = NULL) const = 0;
virtual vec3 eval_hd(const vec3& h, const vec3& d,
const void *user_param = NULL) const;
// evaluate f_r * cos
virtual vec3 evalp(const vec3& i, const vec3& o,
const void *user_param = NULL) const;
virtual vec3 evalp_hd(const vec3& h, const vec3& d,
const void *user_param = NULL) const;
// evaluate f_r * cos / pdf
virtual vec3 evalp_is(float_t u1, float_t u2,
const vec3& o,
vec3 *i, float_t *pdf,
const void *user_param = NULL) const;
// importance sample f_r * cos using two uniform numbers
virtual vec3 sample(float_t u1, float_t u2,
const vec3& o,
const void *user_param = NULL) const;
// evaluate the PDF of a sample
virtual float_t pdf(const vec3& i, const vec3& o,
const void *user_param = NULL) const;
// utilities
static void io_to_hd(const vec3& i, const vec3& o, vec3 *h, vec3 *d);
static void hd_to_io(const vec3& h, const vec3& d, vec3 *i, vec3 *o);
// ctor / dtor
brdf() {}
virtual ~brdf() {}
#if 1 // noncopyable
private:
brdf(const brdf& fr);
brdf& operator=(const brdf& fr);
#endif
};
/* Lambertian BRDF */
class lambert : public brdf {
public:
/* Lambertian Parameters */
class params {
public:
params(const vec3& reflectance = vec3(1));
vec3 m_reflectance;
};
/* Implementation */
vec3 eval(const vec3& i, const vec3& o,
const void *user_param = NULL) const;
};
/* MERL BRDF */
class merl : public brdf {
std::vector<double> m_samples;
public:
merl(const char *filename);
vec3 eval(const vec3& in, const vec3& out,
const void *user_param = NULL) const;
const std::vector<double>& get_samples() const {return m_samples;}
};
/* UTIA BRDF */
class utia : public brdf {
std::vector<double> m_samples;
double m_norm;
public:
utia(const char *filename);
vec3 eval(const vec3& in, const vec3& out,
const void *user_param = NULL) const;
const std::vector<double>& get_samples() const {return m_samples;}
private:
void normalize();
};
/* Fresnel API */
namespace fresnel {
/* Utilities */
void ior_to_f0(float_t ior, float_t *f0);
void ior_to_f0(const vec3& ior, vec3 *f0);
void f0_to_ior(float_t f0, float_t *ior);
void f0_to_ior(const vec3& f0, vec3 *ior);
/* Interface */
class impl {
public:
virtual ~impl() {}
virtual vec3 eval(float_t cos_theta_d) const = 0;
virtual impl *copy() const = 0;
};
/* Ideal Specular Reflection */
class ideal : public impl {
public:
vec3 eval(float_t cos_theta_d) const {return vec3(1);}
impl *copy() const {return new ideal();}
};
/* Fresnel for Unpolarized Light */
class unpolarized : public impl {
vec3 ior; // index of refraction
public:
unpolarized(const vec3& ior): ior(ior) {}
vec3 eval(float_t cos_theta_d) const;
impl *copy() const {return new unpolarized(*this);}
};
/* Schlick's Fresnel */
class schlick : public impl {
vec3 f0; // backscattering fresnel
public:
schlick(const vec3& f0);
vec3 eval(float_t cos_theta_d) const;
impl *copy() const {return new schlick(*this);}
};
/* SGD's Fresnel */
class sgd : public impl {
vec3 f0, f1;
public:
sgd(const vec3& f0, const vec3& f1) : f0(f0), f1(f1) {}
vec3 eval(float_t cos_theta_d) const;
impl *copy() const {return new sgd(*this);}
};
/* Arbitrary Fresnel Function */
class spline : public impl {
std::vector<vec3> m_points;
public:
explicit spline(const std::vector<vec3>& points): m_points(points) {}
const std::vector<vec3>& get_points() const {return m_points;}
vec3 eval(float_t cos_theta_d) const;
impl *copy() const {return new spline(*this);}
};
} // namespace fresnel
/* Microfacet API */
class microfacet : public brdf {
public:
/* microfacet parameters */
class params {
friend class microfacet;
public:
// Factories
static params standard();
static params isotropic(float_t a);
static params elliptic(float_t a1, float_t a2, float_t phi_a = 0.0);
static params pdfparams(float_t ax, float_t ay, float_t rho = 0.0,
float_t tx_n = 0.0, float_t ty_n = 0.0);
// mutators
void set_ellipse(float_t a1, float_t a2, float_t phi_a = 0.0);
void set_pdfparams(float_t ax, float_t ay, float_t rho = 0.0,
float_t tx_n = 0.0, float_t ty_n = 0.0);
void set_location(float_t tx_n, float_t ty_n);
void set_location(const vec3& n);
// accessors
void get_ellipse(float_t *a1, float_t *a2, float_t *phi_a = NULL) const;
void get_pdfparams(float_t *ax, float_t *ay, float_t *rho = NULL,
float_t *tx_n = NULL, float_t *ty_n = NULL) const;
void get_location(float_t *tx_n, float_t *ty_n) const;
void get_location(vec3 *n) const;
// Ctors (prefer factories for construction)
params(float_t a1 = 1.0, float_t a2 = 1.0, float_t phi_a = 0.0);
params(float_t ax, float_t ay, float_t rho, float_t tx_n, float_t ty_n);
private:
vec3 m_n; // mean normal
float_t m_a1, m_a2, m_phi_a; // ellipse parameters
float_t m_ax, m_ay; // scale parameters
float_t m_rho, m_sqrt_one_minus_rho_sqr; // correlation
float_t m_tx_n, m_ty_n; // location parameters
};
// Dtor
virtual ~microfacet() {delete m_fresnel;}
// BRDF interface
vec3 eval(const vec3& i, const vec3& o,
const void *user_param = NULL) const;
vec3 evalp(const vec3& i, const vec3& o,
const void *user_param = NULL) const;
vec3 sample(float_t u1, float_t u2, const vec3& o,
const void *user_param = NULL) const;
float_t pdf(const vec3& i, const vec3& o,
const void *user_param = NULL) const;
vec3 evalp_is(float_t u1, float_t u2, const vec3& o,
vec3 *i, float_t *pdf, const void *user_param = NULL) const;
// eval queries
vec3 fresnel(float_t cos_theta_d) const {return m_fresnel->eval(cos_theta_d);}
float_t ndf(const vec3& h,
const params& params = params::standard()) const;
float_t gaf(const vec3& h, const vec3& i, const vec3& o,
const params& params = params::standard()) const;
float_t g1(const vec3& h, const vec3& k,
const params& params = params::standard()) const;
float_t sigma(const vec3& k,
const params& params = params::standard()) const;
float_t p22(float_t x, float_t y,
const params& params = params::standard()) const;
float_t vp22(float_t x, float_t y, const vec3& k,
const params& params = params::standard()) const;
float_t vndf(const vec3& h, const vec3& k,
const params& params = params::standard()) const;
// sampling queries
virtual bool supports_smith_vndf_sampling() const = 0;
virtual float_t qf2(float_t u, const vec3& k) const;
virtual float_t qf3(float_t u, const vec3& k, float_t qf2) const;
// mutators
void set_shadow(bool shadow) {m_shadow = shadow;}
void set_fresnel(const fresnel::impl& f);
// accessors
int get_shadow() const {return m_shadow;}
const fresnel::impl& get_fresnel() const {return *m_fresnel;}
protected:
// ctor
microfacet(const fresnel::impl& f = fresnel::ideal(),
bool shadow = true);
// microfacet interfaces
virtual float_t sigma_std(const vec3& k) const = 0; // standard masking
virtual float_t p22_std(float_t x, float_t y) const = 0; // standard slope pdf
// sampling functions
virtual void sample_vp22_std_smith(float_t u1, float_t u2, const vec3& k,
float_t *xslope, float_t *yslope) const;
virtual void sample_vp22_std_nmap(float_t u1, float_t u2, const vec3& k,
float_t *xslope, float_t *yslope) const = 0;
// members
const fresnel::impl *m_fresnel;
bool m_shadow;
};
/* Radial microfacets */
class radial : public microfacet {
public:
radial(const fresnel::impl& fresnel = fresnel::ideal(),
bool shadow = true): microfacet(fresnel, shadow)
{}
// queries
virtual float_t p22_radial(float_t r_sqr) const = 0;
virtual float_t sigma_std_radial(float_t cos_theta_k) const = 0;
virtual float_t cdf_radial(float_t r) const = 0;
virtual float_t qf_radial(float_t u) const = 0;
virtual float_t qf2_radial(float_t u,
float_t cos_theta_k,
float_t sin_theta_k) const;
virtual float_t qf3_radial(float_t u, float_t qf2) const;
private:
// eval
float_t p22_std(float_t x, float_t y) const;
float_t sigma_std(const vec3& k) const;
// sample
void sample_vp22_std_nmap(float_t u1, float_t u2, const vec3& o,
float_t *xslope, float_t *yslope) const;
void sample_vp22_std_smith(float_t u1, float_t u2, const vec3& o,
float_t *xslope, float_t *yslope) const;
};
/* Beckmann Microfacet NDF */
class beckmann : public radial {
public:
// Linear Representation
class lrep {
friend class beckmann;
public:
// Ctor
lrep(float_t E1 = 0, float_t E2 = 0,
float_t E3 = 1, float_t E4 = 1,
float_t E5 = 0);
// linear operators
lrep operator+(const lrep& r) const; // combine
lrep operator*(float_t sc) const; // scale
lrep& operator+=(const lrep& rep); // compound combine
lrep& operator*=(float_t sc); // compound scale
// utilities
void scale(float_t sc);
void scale(float_t x, float_t y);
void shear(float_t x, float_t y);
void reparameterize(float_t dudx, float_t dvdx,
float_t dudy, float_t dvdy);
private:
// members
float_t m_E1, m_E2; // first order slope moments
float_t m_E3, m_E4; // second order slope moments
float_t m_E5; // first order joint slope moment
};
// utilities
static void params_to_lrep(const microfacet::params& params, lrep *lrep);
static void lrep_to_params(const lrep& lrep, microfacet::params *params);
// ctor
beckmann(const fresnel::impl& fresnel = fresnel::ideal(),
bool shadow = true): radial(fresnel, shadow)
{}
// queries
float_t p22_radial(float_t r_sqr) const;
float_t sigma_std_radial(float_t cos_theta_k) const;
float_t cdf_radial(float_t r) const;
float_t qf_radial(float_t u) const;
float_t qf1(float_t u) const;
float_t qf2_radial(float_t u,
float_t cos_theta_k, float_t sin_theta_k) const;
float_t qf3_radial(float_t u, float_t qf2) const;
bool supports_smith_vndf_sampling() const {return true;}
};
/* GGX Microfacet NDF */
class ggx : public radial {
public:
ggx(const fresnel::impl& fresnel = fresnel::ideal(),
bool shadow = true): radial(fresnel, shadow)
{}
// queries
float_t p22_radial(float_t r_sqr) const;
float_t sigma_std_radial(float_t cos_theta_k) const;
float_t cdf_radial(float_t r) const;
float_t qf_radial(float_t u) const;
float_t qf1(float_t u) const;
float_t qf2_radial(float_t u,
float_t cos_theta_k, float_t sin_theta_k) const;
float_t qf3_radial(float_t u, float_t qf2) const;
bool supports_smith_vndf_sampling() const {return true;}
private:
float_t qf3_rational_approx(float_t u) const;
};
/* Tabulated Microfacet NDF */
class tabular : public radial {
std::vector<float_t> m_p22;
std::vector<float_t> m_sigma;
// tables for Nmap VNDF sampling
std::vector<float_t> m_cdf;
std::vector<float_t> m_qf;
public:
tabular(const brdf& brdf, int resolution, bool shadow = true);
static microfacet::params fit_beckmann_parameters(const tabular& tabular);
static microfacet::params fit_ggx_parameters(const tabular& tabular);
const std::vector<float_t>& get_p22v() const {return m_p22;}
const std::vector<float_t>& get_sigmav() const {return m_sigma;}
const std::vector<float_t>& get_cdfv() const {return m_cdf;}
const std::vector<float_t>& get_qfv() const {return m_qf;}
// queries
float_t p22_radial(float_t r_sqr) const;
float_t sigma_std_radial(float_t cos_theta_k) const;
float_t cdf_radial(float_t r) const;
float_t qf_radial(float_t u) const;
bool supports_smith_vndf_sampling() const {return false;}
private:
// eval precomputations
void compute_p22_smith(const brdf& brdf, int res);
void compute_fresnel(const brdf& brdf, int res);
void normalize_p22();
void compute_sigma();
// sample precomputations
void compute_cdf();
void compute_qf();
void compute_pdf1();
void normalize_pdf1();
};
/* Tabulated Anisotropic Microfacet NDF */
class tabular_anisotropic : public microfacet {
std::vector<float_t> m_p22;
std::vector<float_t> m_sigma;
std::vector<float_t> m_pdf1;
std::vector<float_t> m_pdf2;
std::vector<float_t> m_cdf1;
std::vector<float_t> m_cdf2;
std::vector<float_t> m_qf1;
std::vector<float_t> m_qf2;
int m_elevation_res;
int m_azimuthal_res;
bool supports_smith_vndf_sampling() const {return false;}
public:
tabular_anisotropic(const brdf& brdf,
int elevation_res,
int azimuthal_res,
bool shadow = true);
static microfacet::params fit_beckmann_parameters(const tabular_anisotropic& tabular);
static microfacet::params fit_ggx_parameters(const tabular_anisotropic& tabular);
const std::vector<float_t>& get_p22v(int *elev_cnt, int *azim_cnt) const;
const std::vector<float_t>& get_sigmav(int *elev_cnt, int *azim_cnt) const;
// queries
float_t pdf1(float_t phi) const;
float_t pdf2(float_t theta, float_t phi) const;
float_t cdf1(float_t phi) const;
float_t cdf2(float_t theta, float_t phi) const;
float_t qf1(float_t u1) const;
float_t qf2(float_t u, float_t phi) const;
private:
// eval interface
float_t sigma_std(const vec3& k) const; // standard masking
float_t p22_std(float_t x, float_t y) const; // standard slope pdf
float_t p22_std_theta_phi(float_t theta, float_t phi) const; // standard slope pdf
// sample interface
void sample_vp22_std_nmap(float_t u1, float_t u2, const vec3& k,
float_t *xslope, float_t *yslope) const;
// eval precomputations
void compute_p22_smith(const brdf& brdf);
void compute_fresnel(const brdf& brdf, int res);
void normalize_p22();
void compute_sigma();
// sample precomputations
void compute_pdf1();
void compute_pdf2();
void normalize_pdf1();
void normalize_pdf2();
void compute_cdf1();
void compute_cdf2();
void compute_qf1();
void compute_qf2();
};
/* Shifted Gamma Distribution BRDF */
class sgd : public brdf {
struct data {
const char *name;
const char *otherName;
double rhoD[3];
double rhoS[3];
double alpha[3];
double p[3];
double f0[3];
double f1[3];
double kap[3];
double lambda[3];
double c[3];
double k[3];
double theta0[3];
double error[3];
};
static const data s_data[100];
const fresnel::impl *m_fresnel;
const data *m_data;
public:
sgd(const char *name);
~sgd() {delete m_fresnel;}
vec3 eval(const vec3& i, const vec3& o,
const void *user_param = NULL) const;
vec3 ndf(const vec3& h) const;
vec3 gaf(const vec3& h, const vec3& i, const vec3& o) const;
vec3 g1(const vec3& k) const;
vec3 fresnel(float_t cos_theta_d) const {return m_fresnel->eval(cos_theta_d);}
const fresnel::impl &get_fresnel() const {return *m_fresnel;}
};
/* ABC Distribution BRDF */
class abc : public brdf {
struct data {
const char *name;
double kD[3];
double A[3];
double B;
double C;
double ior;
};
static const data s_data[100];
const fresnel::impl *m_fresnel;
const data *m_data;
public:
abc(const char *name);
~abc() {delete m_fresnel;}
vec3 eval(const vec3& i, const vec3& o,
const void *user_param = NULL) const;
vec3 ndf(const vec3& h) const;
float_t gaf(const vec3& h, const vec3& i, const vec3& o) const;
vec3 fresnel(float_t cos_theta_d) const {return m_fresnel->eval(cos_theta_d);}
const fresnel::impl &get_fresnel() const {return *m_fresnel;}
};
} // namespace djb
//
//
//// end header file /////////////////////////////////////////////////////
#endif // DJB_INCLUDE_DJ_BRDF_H
#if DJ_BRDF_IMPLEMENTATION
#include <cmath>
#include <cstdarg>
#include <iostream> // std::ios, std::istream, std::cout
#include <fstream> // std::filebuf
#include <cstring> // memcpy
#include <stdint.h> // uint32_t
#ifndef DJB_ASSERT
# include <assert.h>
# define DJB_ASSERT(x) assert(x)
#endif
#ifndef DJB_LOG
# include <stdio.h>
# define DJB_LOG(format, ...) fprintf(stdout, format, ##__VA_ARGS__)
#endif
#ifndef M_PI
# define M_PI 3.1415926535897932384626433832795
#endif
#ifdef _MSC_VER
# pragma warning(disable: 4244) // possible loss of data
#endif
namespace djb {
// *************************************************************************************************
// utility API
template<typename T> static T min(const T& a, const T& b) {return a < b ? a : b;}
template<typename T> static T max(const T& a, const T& b) {return a > b ? a : b;}
template<typename T> static T sat(const T& x) {return min(T(1), max(T(0), x));}
exc::exc(const char *fmt, ...)
{
char buf[256];
va_list args;
va_start(args, fmt);
vsnprintf(buf, 256, fmt, args);
va_end(args);
m_str = std::string(buf);
}
vec3::vec3(float_t theta, float_t phi)
{
float_t s = sin(theta);
x = s * cos(phi);
y = s * sin(phi);
z = cos(theta);
}
#define OP operator
#define V3 vec3
V3 OP*(float_t a, const V3& b) {return V3(a * b.x, a * b.y, a * b.z);}
V3 OP*(const V3& a, float_t b) {return V3(b * a.x, b * a.y, b * a.z);}
V3 OP/(const V3& a, float_t b) {return (1.0 / b) * a;}
V3 OP*(const V3& a, const V3& b) {return V3(a.x * b.x, a.y * b.y, a.z * b.z);}
V3 OP/(const V3& a, const V3& b) {return V3(a.x / b.x, a.y / b.y, a.z / b.z);}
V3 OP+(const V3& a, const V3& b) {return V3(a.x + b.x, a.y + b.y, a.z + b.z);}
V3 OP-(const V3& a, const V3& b) {return V3(a.x - b.x, a.y - b.y, a.z - b.z);}
V3& OP+=(V3& a, const V3& b) {a.x+= b.x; a.y+= b.y; a.z+= b.z; return a;}
V3& OP*=(V3& a, const V3& b) {a.x*= b.x; a.y*= b.y; a.z*= b.z; return a;}
V3& OP*=(V3& a, float_t b) {a.x*= b; a.y*= b; a.z*= b; return a;}
#undef V3
#undef OP
static float_t inversesqrt(float_t x)
{
DJB_ASSERT(x > 0.0);
return (1.0 / sqrt(x));
}
static float_t dot(const vec3& a, const vec3& b)
{
return (a.x * b.x + a.y * b.y + a.z * b.z);
}
static vec3 cross(const vec3& a, const vec3& b)
{
return vec3(a.y * b.z - a.z * b.y,
a.z * b.x - a.x * b.z,
a.x * b.y - a.y * b.x);
}
static vec3 normalize(const vec3& v)
{
float_t mag_sqr = dot(v, v);
if (mag_sqr <= 0.0) printf("%f %f %f\n", v.x, v.y, v.z);
DJB_ASSERT(mag_sqr > 0.0 && "invalid vector magnitude");
return (inversesqrt(mag_sqr) * v);
}
template<typename T>
static T max3(const T& x, const T& y, const T& z)
{
T m = x;
(m < y) && (m = y);
(m < z) && (m = z);
return m;
}
static void xyz_to_theta_phi(const vec3& p, float_t *theta, float_t *phi)
{
if (p.z > 0.99999) {
(*theta) = (*phi) = 0.0;
} else if (p.z < -0.99999) {
(*theta) = M_PI;
(*phi) = 0.0;
} else {
(*theta) = acos(p.z);
(*phi) = atan2(p.y, p.x);
}
}
// *************************************************************************************************
// Special functions
// See: http://www.johndcook.com/blog/cpp_erf/, by John D. Cook
static float_t erf(float_t x)
{
// constants
float_t a1 = (float_t)0.254829592;
float_t a2 = -(float_t)0.284496736;
float_t a3 = (float_t)1.421413741;
float_t a4 = -(float_t)1.453152027;
float_t a5 = (float_t)1.061405429;
float_t p = (float_t)0.3275911;
// Save the sign of x
int sign = 1;
if (x < 0)
sign = -1;
x = fabs(x);
// A&S formula 7.1.26
float_t t = 1.0/(1.0 + p*x);
float_t y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x);
return sign * y;
}
// See: Approximating the erfinv function, by Mike Giles
static float_t erfinv(float_t u)
{
float_t w, p;
w = -logf(((float_t)1.0 - u) * ((float_t)1.0 + u));
if (w < (float_t)5.0) {
w = w - (float_t)2.500000;
p = (float_t)2.81022636e-08;
p = (float_t)3.43273939e-07 + p * w;
p = (float_t)-3.5233877e-06 + p * w;
p = (float_t)-4.39150654e-06 + p * w;
p = (float_t)0.00021858087 + p * w;
p = (float_t)-0.00125372503 + p * w;
p = (float_t)-0.00417768164 + p * w;
p = (float_t)0.246640727 + p * w;
p = (float_t)1.50140941 + p * w;
} else {
w = sqrt(w) - (float_t)3.0;
p = (float_t)-0.000200214257;
p = (float_t)0.000100950558 + p * w;
p = (float_t)0.00134934322 + p * w;
p = (float_t)-0.00367342844 + p * w;
p = (float_t)0.00573950773 + p * w;
p = (float_t)-0.0076224613 + p * w;
p = (float_t)0.00943887047 + p * w;
p = (float_t)1.00167406 + p * w;
p = (float_t)2.83297682 + p * w;
}
return p * u;
}
// *************************************************************************************************
// Sample warps
static void
uniform_to_concentric(float_t u1, float_t u2, float_t *x, float_t *y)
{
/* Concentric map code with less branching (by Dave Cline), see
http://psgraphics.blogspot.ch/2011/01/improved-code-for-concentric-map.html */
float_t r1 = 2.0 * u1 - 1.0;
float_t r2 = 2.0 * u2 - 1.0;
float_t phi, r;
if (r1 == 0 && r2 == 0) {
r = phi = 0;
} else if (r1 * r1 > r2 * r2) {
r = r1;
phi = (M_PI / 4.0) * (r2 / r1);
} else {
r = r2;
phi = (M_PI / 2.0) - (r1 / r2) * (M_PI / 4.0);
}
(*x) = r * cos(phi);
(*y) = r * sin(phi);
}
// *************************************************************************************************
// Geometric operations
//---------------------------------------------------------------------------
// rotate vector along one axis
static vec3 rotate_vector(const vec3& x, const vec3& axis, float_t angle)
{
float_t cos_angle = cos(angle);
float_t sin_angle = sin(angle);
vec3 out = cos_angle * x;
float_t tmp1 = dot(axis, x);
float_t tmp2 = tmp1 * (1.0 - cos_angle);
out+= axis * tmp2;
out+= sin_angle * cross(axis, x);
return out;
}
// *************************************************************************************************
// BRDF API
void brdf::io_to_hd(const vec3& i, const vec3& o, vec3 *h, vec3 *d)
{
const vec3 y_axis = vec3(0, 1, 0);
const vec3 z_axis = vec3(0, 0, 1);
float_t theta_h, phi_h;
(*h) = normalize(i + o);
xyz_to_theta_phi(*h, &theta_h, &phi_h);
vec3 tmp = rotate_vector(i, z_axis, -phi_h);
(*d) = normalize(rotate_vector(tmp, y_axis, -theta_h));
}
void brdf::hd_to_io(const vec3& h, const vec3& d, vec3 *i, vec3 *o)
{
const vec3 y_axis = vec3(0, 1, 0);
const vec3 z_axis = vec3(0, 0, 1);
float_t theta_h, phi_h;
xyz_to_theta_phi(h, &theta_h, &phi_h);
vec3 tmp = rotate_vector(d, y_axis, theta_h);
(*i) = normalize(rotate_vector(tmp, z_axis, phi_h));
(*o) = normalize(2.0 * dot((*i), h) * h - (*i));
}
vec3 brdf::eval_hd(const vec3& h, const vec3& d, const void *user_param) const
{
vec3 i, o;
hd_to_io(h, d, &i, &o);
return eval(i, o, user_param);
}
vec3 brdf::evalp(const vec3& i, const vec3& o, const void *user_param) const
{
return eval(i, o, user_param) * i.z;
}
vec3 brdf::evalp_hd(const vec3& h, const vec3& d, const void *user_param) const
{
vec3 i, o;
hd_to_io(h, d, &i, &o);
return eval(i, o, user_param) * i.z;
}
vec3
brdf::evalp_is(
float_t u1, float_t u2,
const vec3& o, vec3 *i, float_t *pdf,
const void *user_param
) const {
const vec3 i_ = sample(u1, u2, o, user_param);
float_t pdf_ = this->pdf(i_, o);
if (i) (*i) = i_;
if (pdf) (*pdf) = pdf_;
return (evalp(i_, o, user_param) / pdf_);
}
vec3
brdf::sample(float_t u1, float_t u2, const vec3& o, const void *user_param) const
{
float_t x, y, z;
uniform_to_concentric(u1, u2, &x, &y);
z = sqrt(1.0 - x * x - y * y);
DJB_ASSERT(z >= 0.0);
return vec3(x, y, z);
}
float_t brdf::pdf(const vec3& i, const vec3& o, const void *user_param) const
{
return /* cos(theta_i) */i.z / M_PI;
}
// *************************************************************************************************
// Lambertian API implementation
//---------------------------------------------------------------------------
// Lambertian Params
/* Constructors */
lambert::params::params(const vec3& reflectance):
m_reflectance(reflectance)
{}
//---------------------------------------------------------------------------
// Lambertian Evaluation
vec3
lambert::eval(const vec3& i, const vec3& o, const void *user_param) const
{
const lambert::params params =
user_param ? *reinterpret_cast<const lambert::params *>(user_param)
: lambert::params();
return (params.m_reflectance / M_PI);
}
// *************************************************************************************************
// MERL API implementation
// 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
// <[email protected]>.
// 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.
#define MERL_SAMPLING_RES_THETA_H 90
#define MERL_SAMPLING_RES_THETA_D 90
#define MERL_SAMPLING_RES_PHI_D 360
#define MERL_RED_SCALE (1.00 / 1500.0)
#define MERL_GREEN_SCALE (1.15 / 1500.0)
#define MERL_BLUE_SCALE (1.66 / 1500.0)
//---------------------------------------------------------------------------
// Lookup theta_half index
// This is a non-linear mapping!
// In: [0 .. pi/2]
// Out: [0 .. 89]
static int
theta_half_index(float_t theta_half)
{
if (theta_half <= 0.0)
return 0;
float_t theta_half_deg = ((theta_half / (M_PI/2.0)) * MERL_SAMPLING_RES_THETA_H);
float_t temp = theta_half_deg * MERL_SAMPLING_RES_THETA_H;
temp = sqrt(temp);
int ret_val = (int)temp;
if (ret_val < 0)
ret_val = 0;
else if (ret_val >= MERL_SAMPLING_RES_THETA_H)
ret_val = MERL_SAMPLING_RES_THETA_H - 1;
return ret_val;
}
//---------------------------------------------------------------------------
// Lookup theta_diff index
// In: [0 .. pi/2]
// Out: [0 .. 89]
static int
theta_diff_index(float_t theta_diff)
{
int tmp = theta_diff / (M_PI * 0.5) * MERL_SAMPLING_RES_THETA_D;
if (tmp < 0)
return 0;
else if (tmp < MERL_SAMPLING_RES_THETA_D - 1)
return tmp;
else
return MERL_SAMPLING_RES_THETA_D - 1;
}
//---------------------------------------------------------------------------
// Lookup phi_diff index
static int
phi_diff_index(float_t 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 = phi_diff / M_PI * MERL_SAMPLING_RES_PHI_D / 2;
if (tmp < 0)
return 0;
else if (tmp < MERL_SAMPLING_RES_PHI_D / 2 - 1)
return tmp;
else
return MERL_SAMPLING_RES_PHI_D / 2 - 1;
}
// XXX End of
// Copyright 2005 Mitsubishi Electric Research Laboratories All Rights Reserved.
//---------------------------------------------------------------------------
// MERL Contructor
merl::merl(const char *filename)
{
std::fstream f(filename, std::fstream::in | std::fstream::binary);
int32_t n, dims[3];
// check file
if (!f.is_open())
throw exc("djb_error: Failed to open %s\n", filename);
// read header
f.read((char *)dims, /*bytes*/4 * 3);
n = dims[0] * dims[1] * dims[2];
if (n <= 0)
throw exc("djb_error: Failed to read MERL header\n");
// allocate brdf and read data
m_samples.resize(3 * n);
f.read((char *)&m_samples[0], sizeof(double) * 3 * n);
if (f.fail())
throw exc("djb_error: Reading %s failed\n", filename);
}
//---------------------------------------------------------------------------
// look up the BRDF.
vec3 merl::eval(const vec3& i, const vec3& o, const void *user_param) const
{
// convert to half / diff angle coordinates
vec3 h, d;
float_t theta_h, phi_h, theta_d, phi_d;
io_to_hd(i, o, &h, &d);
xyz_to_theta_phi(h, &theta_h, &phi_h);
xyz_to_theta_phi(d, &theta_d, &phi_d);
// compute indexes
int idx_r = phi_diff_index(phi_d)
+ theta_diff_index(theta_d)
* MERL_SAMPLING_RES_PHI_D / 2
+ theta_half_index(theta_h)