HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
UT_Cdf.h
Go to the documentation of this file.
1 /*
2  * PROPRIETARY INFORMATION. This software is proprietary to
3  * Side Effects Software Inc., and is not to be reproduced,
4  * transmitted, or disclosed in any way without written permission.
5  *
6  * NAME: UT_Cdf.h ( UT Library, C++)
7  *
8  * COMMENTS: An N-Dimensional cumulative distribution function (CDF)
9  * for sampling.
10  */
11 
12 #ifndef __UT_Cdf__
13 #define __UT_Cdf__
14 
15 #include "UT_API.h"
16 #include <SYS/SYS_Types.h>
17 #include <SYS/SYS_Math.h>
18 #include "UT_Assert.h"
19 #include "UT_Algorithm.h"
20 #include "UT_RootFinder.h"
21 #include "UT_VoxelArray.h"
22 
24 {
25  template<class InputIt>
26  static inline int eval(InputIt cdf_begin,
27  const int &res,
28  const float &u)
29  {
30  int i = 0;
31  for ( ; i < res; ++i)
32  {
33  if (cdf_begin[i] > u) { break; }
34  }
35 
36  return i;
37  }
38 };
39 
41 {
42  template<class InputIt>
43  static inline int eval(InputIt cdf_begin,
44  const int &res,
45  const float &u)
46  {
47  InputIt u_i = std::upper_bound(cdf_begin, cdf_begin+res, u);
48 
49  return SYSmax(int(0), int(u_i-cdf_begin));
50  }
51 };
52 
54 {
55  template<class InputIt>
56  static inline int eval(InputIt cdf_begin,
57  const int &res,
58  const float &u)
59  {
60  InputIt u_i = UTfastUpperBound(cdf_begin, cdf_begin+res, u);
61 
62  return SYSmax(int(0), int(u_i-cdf_begin));
63  }
64 };
65 
67 
68 template<int N>
70 {
71  static inline int eval(const int *res,
72  const int *i)
73  {
74  return 0;
75  }
76 };
77 
78 template<>
79 struct UT_CdfIndexer<1>
80 {
81  static inline int eval(const int *res,
82  const int *i)
83  {
84  return i[0];
85  }
86 };
87 
88 template<>
89 struct UT_CdfIndexer<2>
90 {
91  static inline int eval(const int *res,
92  const int *i)
93  {
94  return i[1] * res[0] + i[0];
95  }
96 };
97 
98 template<>
99 struct UT_CdfIndexer<3>
100 {
101  static inline int eval(const int *res,
102  const int *i)
103  {
104  return (i[0] * res[1] + i[1]) * res[2] + i[2];
105  }
106 };
107 
109 {
110  enum { N = 1 };
111  inline float operator()(const int *i) const { return 0.0F; }
112 };
113 
114 template<int N>
116 {
117  template<class F>
118  static inline void eval(float *cdf,
119  const int *res,
120  const F &ftor)
121  {
122  }
123 };
124 
125 template<>
127 {
128  enum { N = 1 };
129 
130  template<class F>
131  static inline void eval(float *cdf,
132  const int *res,
133  const F &ftor)
134  {
135  int i[N] = {0};
136  int j = 0;
137  for ( ; i[0] < res[0]; ++i[0])
138  {
139  cdf[j++] = ftor(i);
140  }
141  }
142 };
143 
144 template<>
146 {
147  enum { N = 2 };
148 
149  template<class F>
150  static inline void eval(float *cdf,
151  const int *res,
152  const F &ftor)
153  {
154  int i[N] = {0, 0};
155  int j = 0;
156  for (i[0] = 0; i[0] < res[0]; ++i[0])
157  {
158  for (i[1] = 0; i[1] < res[1]; ++i[1])
159  {
160  cdf[j++] = ftor(i);
161  }
162  }
163  }
164 };
165 
166 template<>
168 {
169  enum { N = 3 };
170 
171  template<class F>
172  static inline void eval(float *cdf,
173  const int *res,
174  const F &ftor)
175  {
176  int i[N] = {0, 0, 0};
177  int j = 0;
178  for (i[0] = 0; i[0] < res[0]; ++i[0])
179  {
180  for (i[1] = 0; i[1] < res[1]; ++i[1])
181  {
182  for (i[2] = 0; i[2] < res[2]; ++i[2])
183  {
184  cdf[j++] = ftor(i);
185  }
186  }
187  }
188  }
189 };
190 
192 {
193  enum { N = 3 };
194 
196  : myVoxels(voxels) {}
197 
198  inline float operator()(const int *i) const
199  {
200  return myVoxels->operator()(i[0], i[1], i[2]);
201  }
202 
204 };
205 
206 /// Fill a linear cdf with res entries of a linear gradient.
207 template<class InputIt>
208 inline void
209 UTfillCdfIdentity(InputIt cdf, const int &res)
210 {
211  if (res == 0) { return; }
212  float s = 1.0F / res;
213  for (int i = 0; i < res; ++i)
214  {
215  *cdf++ = s * i;
216  }
217 }
218 
219 /// Sample from a linear cdf with res entries. The return value is the
220 /// sampled array index.
221 template<class S, class InputIt>
222 inline int
223 UTsampleCdf(InputIt cdf_begin, const int &res, const float &u)
224 {
225  return SYSclamp(S::eval(cdf_begin, res, u), 0, res-1);
226 }
227 
228 /// Sample from a linear cdf with res entries. The return value is the
229 /// sampled array index and dval is the offset.
230 template<class S, class InputIt>
231 inline int
232 UTsampleCdf(InputIt cdf_begin, const int &res, const float &u,
233  float &dval, float &pdf)
234 {
235  fpreal range;
236  int idx;
237 
238  dval = u;
239  dval *= cdf_begin[res-1];
240  idx = UTsampleCdf<S>(cdf_begin, res, dval);
241 
242  range = cdf_begin[idx];
243  pdf = 0;
244  if (idx > 0)
245  {
246  dval -= cdf_begin[idx-1];
247  range -= cdf_begin[idx-1];
248  }
249  if (range > 0)
250  {
251  dval /= range;
252  pdf = res*range / cdf_begin[res-1];
253  }
254 
255  return idx;
256 }
257 
258 /// Sample the pdf from a linear cdf at index i with res entries.
259 template<class InputIt>
260 inline void
261 UTsampleCdf(InputIt cdf_begin, const int &res, const int &i, float &pdf)
262 {
263  fpreal range;
264 
265  range = cdf_begin[i];
266  pdf = 0;
267  if (i > 0)
268  {
269  range -= cdf_begin[i-1];
270  }
271  if (range > 0)
272  {
273  pdf = res*range / cdf_begin[res-1];
274  }
275 }
276 
277 /// Find the value x from 0 to 1 where the area from 0 to x under the
278 /// trapezoid (0,0),(1,0),(1,d1),(0,d0) is portion u of the total area.
279 inline float
280 UTsampleTrapezoidPdf(float u, float d0, float d1)
281 {
282  // Solve for x:
283  // 0 = (d1-d0)*x^2 + 2*d0*x - (d0+d1)*u
284  float x0; float x1;
285  int nroots = UT_RootFinder::quadratic(d1-d0, 2*d0, -(d0+d1)*u, x0, x1);
286  // If there are multiple roots, pick the one that is in
287  // or just outside of the [0,1] interval.
288  float x = (nroots == 1 || SYSabs(x0-0.5f) <= SYSabs(x1-0.5f)) ? x0 : x1;
289  return x;
290 }
291 
292 /// Create a PDF from a range of values.
293 template<class InputIt, class OutputIt>
294 inline void
295 UTcreatePdf(InputIt values_begin,
296  InputIt values_end,
297  OutputIt pdf_begin)
298 {
299  typedef typename std::iterator_traits<OutputIt>::value_type output_t;
300 
301  UTnormalizeArray(values_begin, values_end, pdf_begin);
302 
303  output_t sum = 0;
304  std::for_each(pdf_begin, pdf_begin + (values_end-values_begin),
305  [&sum](const output_t &v)
306  {
307  sum += v;
308  });
309 
310  float area = 1.0F / sum;
311  std::transform(pdf_begin, pdf_begin + (values_end-values_begin), pdf_begin,
312  [&area](const output_t &v) { return v * area; });
313 }
314 
315 /// Create a CDF from a range of values.
316 template<class InputIt, class OutputIt>
317 inline void
318 UTcreateCdf(InputIt values_begin,
319  InputIt values_end,
320  OutputIt cdf_begin)
321 {
322  int i = 0;
323  while (values_begin != values_end)
324  {
325  *cdf_begin = *values_begin;
326  if (i > 0) { *cdf_begin += *(cdf_begin-1); }
327  ++values_begin;
328  ++cdf_begin;
329  ++i;
330  }
331 }
332 
333 template<int N = 1,
334  class S = UT_Cdf_default_search_t>
336 {
337  static_assert((N > 0 && N < 4), "unsupported number of dimensions");
338 
339 public:
340  enum { DIM = N };
341 
342  inline UT_Cdf();
343 
344  /// Construct from an array
345  /// data = row-major pdf/cdf data
346  /// res = resolution for each dimension
347  /// data values less than zero will be treated as zero.
348  inline UT_Cdf(float *data,
349  const int *res,
350  const bool &copydata);
351 
352  /// Construct from a value functor
353  /// ftor = value functor
354  /// res = resolution for each dimension
355  /// data values less than zero will be treated as zero.
356  template<class F>
357  inline UT_Cdf(const F &ftor,
358  const int *res);
359 
360  /// Create a cdf without data. You should then use getRawArray()
361  /// or setRawArray() functions and build() to construct the cdf.
362  /// This method avoids duplication of memory for the pdf and cdf
363  /// when only the cdf is needed.
364  /// res = resolution for each dimension
365  inline UT_Cdf(const int *res);
366 
367  inline ~UT_Cdf();
368 
369  /// Clear associated storage.
370  inline void clear();
371 
372  /// Resize storage according to resolution.
373  /// res = resolution for each dimension
374  inline void resize(const int *res);
375 
376  /// Set the data array from an array.
377  /// data = row-major pdf/cdf data
378  /// copydata = should the data be copied?
379  inline void setRawArray(float *data,
380  const bool &copydata);
381 
382  /// Set the data array from a value functor.
383  /// ftor = value functor
384  template<class F>
385  inline void setRawArray(const F &ftor);
386 
387  /// Retrieve the raw data array in row-major order. There will be
388  /// product(res) entries in the array to fill out.
389  inline float *getRawArray() { return myCdf[N-1]; }
390 
391  /// Build the cdf from the row-major data array. If you provided data
392  /// to the constructor, you don't need to call this method unless
393  /// you've modified the data through getRawArray().
394  inline void build();
395 
396  /// Draw an importance sample from the cdf. u and dval must have
397  /// 'N' entries. In this version, dval is the resulting sample value
398  /// for each dimension in the range [0, 1).
399  inline void sample(const float *u, float *dval) const;
400  inline void sample(const float *u, float *dval, float &pdf) const;
401 
402  /// Same as above but produces an exact index.
403  inline void sample(const float *u, int *didx) const;
404 
405  /// Same as above but produces an exact index. (dval) above relates to
406  /// (didx, doff) using the following identity:
407  /// dval[dim] = (didx[dim] + doff[dim]) / res[dim]
408  inline void sample(const float *u, int *didx, float *doff, float &pdf) const;
409 
410  /// Evaluate the pdf at the given index.
411  inline void evaluatePdf(const float *u, float &pdf) const;
412 
413  /// Get the sum over all image values
414  inline const float &getSum() const { return myCdf[0][myRes[0]-1]; }
415  inline const float &getISum() const { return myISum; }
416 
417  /// Get the average image value
418  inline const float &getScaling() const { return myScale; }
419  inline const float &getIScaling() const { return myIScale; }
420 
421  inline int getDim() const { return N; }
422  inline const int *getRes() const { return myRes; }
423  inline const int &getRes(int dim) const { return myRes[dim]; }
424 
425  inline void dump() const;
426 
427 private:
428  inline void init();
429  inline void realloc(const int *res);
430 
431  float *myCdf[N];
432  int myRes[N];
433  float myIRes[N];
434  float myISum;
435  float myScale;
436  float myIScale;
437  int myCapacity;
438  int mySize;
439  bool myOwnData;
440 };
441 
442 template<int N, class S>
443 inline
445 {
446  init();
447 }
448 
449 template<int N, class S>
450 inline
452  const int *res,
453  const bool &copydata)
454 {
455  init();
456  myOwnData = copydata;
457  resize(res);
458  setRawArray(data, copydata);
459  build();
460 }
461 
462 template<int N, class S>
463 template<class F>
464 inline
465 UT_Cdf<N, S>::UT_Cdf(const F &ftor,
466  const int *res)
467 {
468  init();
469  resize(res);
470  setRawArray(ftor);
471  build();
472 }
473 
474 template<int N, class S>
475 inline
476 UT_Cdf<N, S>::UT_Cdf(const int *res)
477 {
478  init();
479  resize(res);
480 }
481 
482 template<int N, class S>
483 inline
485 {
486  clear();
487 }
488 
489 template<int N, class S>
490 inline void
492 {
493  for (int i = 0; i < N; ++i)
494  {
495  myCdf[i] = nullptr;
496  myRes[i] = 0;
497  myIRes[i] = 0;
498  }
499  myISum = 0;
500  myScale = 0;
501  myIScale = 0;
502  myCapacity = 0;
503  mySize = 0;
504  myOwnData = true;
505 }
506 
507 template<int N, class S>
508 inline void
510 {
511  if (!myOwnData) { return; }
512  for (int i = 0; i < N; ++i)
513  {
514  if (myCdf[i])
515  {
516  delete[] myCdf[i];
517  myCdf[i] = nullptr;
518  }
519  }
520  myCapacity = 0;
521 }
522 
523 template<int N, class S>
524 inline void
525 UT_Cdf<N, S>::realloc(const int *res)
526 {
527  if (!myOwnData) { return; }
528  clear();
529  myCapacity = 1;
530  for (int i = 0; i < N; ++i)
531  {
532  myCapacity *= res[i];
533  myCdf[i] = new float[myCapacity];
534  }
535 }
536 
537 template<int N, class S>
538 inline void
539 UT_Cdf<N, S>::resize(const int *res)
540 {
541  mySize = 1;
542  for (int i = 0; i < N; ++i)
543  {
544  myRes[i] = res[i];
545  myIRes[i] = 1.0F / myRes[i];
546  mySize *= myRes[i];
547  }
548  if (mySize > myCapacity) { realloc(res); }
549 }
550 
551 template<int N, class S>
552 inline void
554  const bool &copydata)
555 {
556  if (myOwnData && !copydata)
557  {
558  clear();
559  }
560  myOwnData = copydata;
561  if (copydata)
562  {
563  memcpy(myCdf[N-1], data, mySize*sizeof(float));
564  }
565  else
566  {
567  myCdf[N-1] = data;
568  }
569 }
570 
571 template<int N, class S>
572 template<class F>
573 inline void
575 {
576  static_assert(std::is_class<F>::value, "invalid functor object");
577  static_assert(F::N == N, "unsupported number of dimensions");
578 
579  myOwnData = true;
580  UT_CdfValueFill<N>::eval(myCdf[N-1], myRes, ftor);
581 }
582 
583 template<int N, class S>
584 inline void
586 {
587  // Compute the sums for efficient sampling
588  int size = mySize;
589  for (int i = N; i-- > 0; )
590  {
591  size /= myRes[i];
592 
593  int idx = 0;
594  for (int j = 0; j < size; ++j)
595  {
596  // We need non-negative, so we clamp negative to zero.
597  myCdf[i][idx] = SYSmax(myCdf[i][idx], 0.0F);
598  ++idx;
599  for (int k = 1; k < myRes[i]; ++k, ++idx)
600  {
601  myCdf[i][idx] = SYSmax(myCdf[i][idx], 0.0F);
602  myCdf[i][idx] += SYSmax(myCdf[i][idx-1], 0.0F);
603  }
604 
605  if (i > 0)
606  myCdf[i-1][j] = myCdf[i][idx-1];
607  }
608  }
609 
610  UT_ASSERT(size == 1);
611 
612  myScale = getSum();
613  for (int i = 0; i < N; ++i)
614  myScale /= (float)myRes[i];
615 
616  myISum = 1.0F / getSum();
617  myIScale = 1.0F / getScaling();
618 }
619 
620 template<int N, class S>
621 inline void
622 UT_Cdf<N, S>::sample(const float *u, float *dval) const
623 {
624  float pdf = 0;
625 
626  sample(u, dval, pdf);
627 }
628 
629 template<int N, class S>
630 inline void
631 UT_Cdf<N, S>::sample(const float *u, float *dval, float &pdf) const
632 {
633  int didx[N];
634 
635  sample(u, didx, dval, pdf);
636  for (int i = 0; i < N; ++i)
637  {
638  dval[i] += didx[i];
639  dval[i] *= myIRes[i];
640  }
641 }
642 
643 template<int N, class S>
644 inline void
645 UT_Cdf<N, S>::sample(const float *u, int *didx) const
646 {
647  const float *cdf = nullptr;
648  int off = 0;
649  int idx = 0;
650  for (int i = 0; i < N; ++i)
651  {
652  off += idx;
653  off *= myRes[i];
654  cdf = myCdf[i] + off;
655 
656  idx = didx[i] = UTsampleCdf<S>(
657  cdf, myRes[i], u[i] * cdf[myRes[i]-1]);
658  }
659 }
660 
661 template<int N, class S>
662 inline void
663 UT_Cdf<N, S>::sample(const float *u, int *didx, float *doff, float &pdf) const
664 {
665  const float *cdf = nullptr;
666  int off = 0;
667  int idx = 0;
668  float onepdf = 0;
669  pdf = 1;
670  for (int i = 0; i < N; ++i)
671  {
672  off += idx;
673  off *= myRes[i];
674  cdf = myCdf[i] + off;
675 
676  idx = didx[i] = UTsampleCdf<S>(
677  cdf, myRes[i], u[i], doff[i], onepdf);
678  pdf *= onepdf;
679  }
680 }
681 
682 template<int N, class S>
683 inline void
684 UT_Cdf<N, S>::evaluatePdf(const float *u, float &pdf) const
685 {
686  int off = 0;
687  int idx = 0;
688  for (int i = 0; i < N; ++i)
689  {
690  idx = (int)(u[i]*myRes[i]);
691  idx = SYSclamp(idx, 0, myRes[i]-1);
692 
693  off *= myRes[i];
694  off += idx;
695  }
696 
697  pdf = myCdf[N-1][off];
698  if (off > 0)
699  {
700  pdf -= myCdf[N-1][off-1];
701  }
702 
703  // Divide by the average value to normalize the pdf to 1
704  pdf *= getIScaling();
705 }
706 
707 template<int N, class S>
708 inline void
710 {
711  for (int dim = 0; dim < N; ++dim)
712  {
713  fprintf(stderr, "dim %d:\n", dim);
714  for (int i = 0; i < myRes[dim]; ++i)
715  fprintf(stderr, "%d: %f\n", i, myCdf[dim][i]);
716  }
717 }
718 
719 #endif
static void eval(float *cdf, const int *res, const F &ftor)
Definition: UT_Cdf.h:118
void sample(const float *u, float *dval) const
Definition: UT_Cdf.h:622
#define SYSmax(a, b)
Definition: SYS_Math.h:1365
const int & getRes(int dim) const
Definition: UT_Cdf.h:423
static int eval(const int *res, const int *i)
Definition: UT_Cdf.h:71
GLenum GLint * range
Definition: glcorearb.h:1924
static void eval(float *cdf, const int *res, const F &ftor)
Definition: UT_Cdf.h:131
float * getRawArray()
Definition: UT_Cdf.h:389
const float & getSum() const
Get the sum over all image values.
Definition: UT_Cdf.h:414
const GLdouble * v
Definition: glcorearb.h:836
#define SYSabs(a)
Definition: SYS_Math.h:1367
static int eval(const int *res, const int *i)
Definition: UT_Cdf.h:91
float operator()(const int *i) const
Definition: UT_Cdf.h:198
int UTsampleCdf(InputIt cdf_begin, const int &res, const float &u)
Definition: UT_Cdf.h:223
png_uint_32 i
Definition: png.h:2877
uint64 value_type
Definition: GA_PrimCompat.h:29
void UTfillCdfIdentity(InputIt cdf, const int &res)
Fill a linear cdf with res entries of a linear gradient.
Definition: UT_Cdf.h:209
GLsizeiptr size
Definition: glcorearb.h:663
UT_CdfStdBinarySearch UT_Cdf_default_search_t
Definition: UT_Cdf.h:66
static int eval(const int *res, const int *i)
Definition: UT_Cdf.h:81
const float & getIScaling() const
Definition: UT_Cdf.h:419
const float & getISum() const
Definition: UT_Cdf.h:415
GLfloat f
Definition: glcorearb.h:1925
static void eval(float *cdf, const int *res, const F &ftor)
Definition: UT_Cdf.h:172
void dump() const
Definition: UT_Cdf.h:709
#define UT_ASSERT(ZZ)
Definition: UT_Assert.h:102
float UTsampleTrapezoidPdf(float u, float d0, float d1)
Definition: UT_Cdf.h:280
const float & getScaling() const
Get the average image value.
Definition: UT_Cdf.h:418
void setRawArray(float *data, const bool &copydata)
Definition: UT_Cdf.h:553
UT_CdfValueVoxelArrayF(const UT_VoxelArrayF *voxels)
Definition: UT_Cdf.h:195
void evaluatePdf(const float *u, float &pdf) const
Evaluate the pdf at the given index.
Definition: UT_Cdf.h:684
~UT_Cdf()
Definition: UT_Cdf.h:484
static int eval(const int *res, const int *i)
Definition: UT_Cdf.h:101
GLboolean * data
Definition: glcorearb.h:130
const int * getRes() const
Definition: UT_Cdf.h:422
GA_API const UT_StringHolder transform
static int eval(InputIt cdf_begin, const int &res, const float &u)
Definition: UT_Cdf.h:26
void UTcreateCdf(InputIt values_begin, InputIt values_end, OutputIt cdf_begin)
Create a CDF from a range of values.
Definition: UT_Cdf.h:318
void UTcreatePdf(InputIt values_begin, InputIt values_end, OutputIt pdf_begin)
Create a PDF from a range of values.
Definition: UT_Cdf.h:295
UT_Cdf()
Definition: UT_Cdf.h:444
static int eval(InputIt cdf_begin, const int &res, const float &u)
Definition: UT_Cdf.h:43
int getDim() const
Definition: UT_Cdf.h:421
GLsizei const GLfloat * value
Definition: glcorearb.h:823
double fpreal
Definition: SYS_Types.h:269
void UTnormalizeArray(InputIt begin, InputIt end, OutputIt d_begin)
Normalize an array.
Definition: UT_Algorithm.h:282
typedef int
Definition: png.h:1175
void build()
Definition: UT_Cdf.h:585
#define UT_API_TMPL
Definition: UT_API.h:13
GLint GLenum GLint x
Definition: glcorearb.h:408
const UT_VoxelArrayF * myVoxels
Definition: UT_Cdf.h:203
GA_API const UT_StringHolder N
void clear()
Clear associated storage.
Definition: UT_Cdf.h:509
Definition: UT_Cdf.h:335
static void eval(float *cdf, const int *res, const F &ftor)
Definition: UT_Cdf.h:150
static int quadratic(T a, T b, T c, T &v0, T &v1)
void resize(const int *res)
Definition: UT_Cdf.h:539
ForwardIt UTfastUpperBound(ForwardIt first, ForwardIt last, const T &value)
Definition: UT_Algorithm.h:324
float operator()(const int *i) const
Definition: UT_Cdf.h:111
GA_API const UT_StringHolder area
static int eval(InputIt cdf_begin, const int &res, const float &u)
Definition: UT_Cdf.h:56