30:Tf36,md`# Sampling Sampling, the process of creating samples, is often used in rendering. For example, in Ray tracing we sample a point on the pixel for generating a ray. Most of us probably chose a fixed point, the center of the pixel. If you remembered the ray computation we used coordinates (0.5,0.5) to get to the center. Instead of one sample if we wanted four samples we will divide pixel into four equal area sub-pixels and pixel square area into four smaller square areas and take the center of the sub-pixels as the sample coordinates. That means we will use (0.25,0.25), (0.25,0.75), (0.75,0.25) and (0.75,0.75) as the coordinates. So we used a fixed sampling approach, and particularly we used fixed 2D sampling. ### Random Sampling: Instead, we may choose random sampling approach. That means instead of choosing a fixed value for each dimension in the range [0..1] we will randomly choose in that same range. To get such a random value, we will use a random number generator. Such a random number generator is believed to generate numbers uniformly distributed between [0..1]. Unfortunately, the uniform distribution property is valid when the number of samples generated are very large. In practice we use tend only a very small number of samples, say 2 or 4 or 8. When samples are drawn in such a small number, the samples are almost never uniformly distributed within the range. It is not uncommon to find the samples clumped in a smaller region of the range. To avoid such a problem, a *stratified sampling* approach is used, in which the whole range is equally divided into number of required samples, and within the subrange a sample is chosen completel in random, is randomly jittered from the center of the subrange. ### Random Number Generator All random samplings rely on a uniform random number generator. Though most CPU computing languages provide such a generator as a part of their standard tools, they may be missing in a GPU based computing language. So we must write our own generator. Many random number generating algorithms exist in the literature ([see Wiki](https://en.wikipedia.org/wiki/Random_number_generation#Computational_methods)). Here is a uniform random sampler algorithm from Nokia Raytracer (from their WebCL based demo. Unfortunately the link to the page is no more available): ~~~ float random(unsigned int *seed0, unsigned int *seed1) { *seed0 = 36969 * ((*seed0) & 65535) + ((*seed0) >> 16); *seed1 = 18000 * ((*seed1) & 65535) + ((*seed1) >> 16); unsigned int ires = ((*seed0) << 16) + (*seed1); //Convert to float union { float f; unsigned int ui; } res; res.ui = (ires & 0x007fffff) | 0x40000000; return (res.f - 2.) / 2.; } ~~~ The function takes two integer seeds as input and returns a floating point random number in the range 0 to 1. The function also changes the seed values during the computation and the changed seeds must be retained for subsequent use. That is why the seed parameters are passed by reference. We may use this or any similar function for sampling needs in GPU. The seeds must be chosen in random to start with. We will generate the seeds in CPU using an existing random number generator. For example for every pixel we may predefine the seeds in CPU as follows and send it to GPU (an array for all pixels) along with other data: ~~~ maxint = 2147483647; seed0 = Math.random() * maxint; if(seed0 < 2) seed0 = 2; seed1 = Math.random() * maxint; if(seed1 < 2) seed1 = 2; ~~~ The seed array will be updated to GPU during every invocation of the kernel/shader program. So seed data has to be loaded to GPU only once and used in any number of invocation of the kernel/shader. The random number generator for [ShaderToy](https://www.shadertoy.com/view/tsf3Dn) generated the seed for each pixel from the pixel coordinatte, and uses only one seed! You may try the *nextFloat* function from that page. `31:T62b,md`## Sampling a unit line or a unit rectangle By unit line I mean, 0 ≤ x ≤ 1 and unit rectangle I mean 0 ≤ x,y ≤ 1. One may consider *x* to be scalar parameters that define a straight line segment , and *x* and *y* to be scalar parameters that define a rectangle or even parallelelogram. The following examples describe sampling a known number of samples in 1D the range [0..1] and in 2D in the range [0..1] x [0..1]. ###Fixed Sampling in 1D ~~~ for (int i = 0; i < nSamples; ++i) { sample[i] = (i+0.5)/nSamples; } ~~~ ### Stratified Sampling in 1D ~~~ for (int i = 0; i < nSamples; ++i{ sample[i] = min((i + random(&seed0, &seed1))/nSamples); } ~~~ Note that for these two sampling methods one must know a priori exactly how many samples (i.e. *nSamples*) are required. ### Random Sampling in 1D ~~~ for (int i = 0; i < nSamples; ++i) { sample[i] = random(&seed0, &seed1); } ~~~ This latter sampling method may be used for any number of times, without knowing exactly how many samples one wishes to draw. ### Stratified Sampling in 2D ~~~ for (int y = 0; y < nSamples_y; ++y) for (int x = 0; x < Samples_x; ++x) { sample[i].x = (x + random(&seed0, &seed1))/nSamples_x; sample[i].y = (y + random(&seed0, &seed1))/nSamples_y; } } ~~~ ### Random Sampling in 2D ~~~ for (int y = 0; y < nSamples_y; ++y) for (int x = 0; x < Samples_x; ++x) { sample[i].x = random(&seed0, &seed1); sample[i].y = random(&seed0, &seed1); } } ~~~ Like in the 1D case this latter sampling method may be used without prior knowledge of how many samples are required. `32:T6ae,md`### Sampling of Pixel for Ray Generation. We have so far traced a single ray through the center of the pixel. So we used the following equation to compute the ray direction through pixel (i,j): ~~~ ray.direction = -w +(-0.5+(i+0.5)/Cols))*windowWidth*U + (-0.5 +(j+0.5)/Rows)*windowHeight*V ~~~ In the equation we added offset 0.5 to both i and j (i.e. we used(i+0.5) and (j+0.5)) to get to the center of the pixel (i,j). If we wish to send a ray a random location inside the pixel instead, we will simply substitute the fixed offset 0.5 by a random offset to i and a random offset to j. The random offsets must be uniformly sampled between 0 and 1. ~~~ float ξ_1 = random(&seed0, &seed1); float ξ_2 = random(&seed0, &seed1); ray.direction = -w +(-0.5+(i+ξ_1)/Cols))*windowWidth*U + (-0.5 +(j+ξ_2)/Rows)*windowHeight*V ~~~ We will repeat this method to compute any number of random samples through the pixel. Note that we did not use stratified sampling here. If the number of samples are small then samples may not get uniformly distributed inside the pixel area, say for example, samples may get crowded into one corner of the pixel. If we know the number of sampled *a priori* then stratified sampling will provide a little better distribution of samples inside the area. For example, if we wish to send MxN uniformly distributed rays through each pixel then the direction of the ray sample (m,n) though pixel(i,j) will be ~~~ float ξ_1 = (m+random(&seed0, &seed1))/M; float ξ_2 = (n+random(&seed0, &seed1))/N; ray.direction = -w +(-0.5+(i+ξ_1)/Cols))*windowWidth*U + (-0.5 +(j+ξ_2)/Rows)*windowHeight*V ~~~ where *m* takes value in the range (0 .. M-1) and n in the rannge (0 .. N-1). `33:T7f5,md`### Sampling Disk Disk sampling has a number of applications in ray tracing. For example: sampling a thin lens, sampling a disc shaped light source, uniform sampling a hemisphere, etc.. We describe here the three sampling methods for disc sampling. These methods may be used to generate random samples without a priori knowledge of how many samples are needed. For these methods we will assume that a disk of radius 1 is centered around the origin on XY plane. So the sampling computes (x,y) coordinates for such a simple disk. To generated samples on an arbitrary, disk of radius R lying on UV plane and centered around a point C the samples must be transformed using the following equation: ~~~ samplePoint = C + x*R*U + y*R*V ~~~ where U and V are normalized vectors that define the axes of the 2D plane of the disk. #### Rejection sampling ~~~ do { x = 2. * random(&seed0,&seed1) - 1.; // Converting [0..1] range to [-1..1] y = 2. * random(&seed0,&seed1) - 1.; } while (x*x + y*y > 1.); ~~~ #### Uniform area sampling The following method is based on polar coordinates to Cartesian coordinate conversion. It computes θ in the range (0 to 2*π). ~~~ float r = sqrt(random(&seed0,&seed1)); float θ = 2.0 * π * random(&seed0,&seed1); x = r * cos(θ); y = r * sin(θ); ~~~ #### Concentric Disc sampling The following method computes theta in the range (0 to 2*π). ~~~ // Map uniform random number in [0..1] range to [-1,1] float sx = 2 * random(&seed0,&seed1) - 1; float sy = 2 * random(&seed0,&seed1) - 1; // Map square to (r,θ) if (sx == 0.0 && sy == 0.0) { x = 0.0; y = 0.0; } else if (sx >= -sy) { if (sx > sy) { // Handle first region of disk r = sx; if (sy > 0.0) θ = sy/r; else θ = 8.0 + sy/r; } else { // Handle second region of disk r = sy; θ = 2.0 - sx/r; } } else { if (sx <= sy) { // Handle third region of disk r = -sx; θ = 4.0 - sy/r; } else { // Handle fourth region of disk r = -sy; θ = 6.0 + sx/r; } } θ *= π / 4.; x = r * cos(θ); y = r * sin(θ); ~~~ `34:T66a,md`### Uniformly Sampling of Sphere There are many ways to uniformly sample the sphere. We describe a method based on the following principles: The area of a sphere equals the area of every right circular cylinder circumscribed about the sphere excluding the bases. One can verify this simple principle. ~~~ Area of the sphere is 4*π*R*R. ~~~ The height of the right circular cylinder circumscribed about the sphere is 2R. ~~~ The area of the cylinder of height 2*R (without the base and top) is 2*π*R*height = 2*π*R*2*R = 4*π*R*R. ~~~ So if we uniformly sample the area of the right circular cylinder described above then we will get a uniform sampling of the sphere. Uniformly sampling the cylinder surface is simple. So algorithm is - Uniformly sample the height between +R to -R. - Uniformly sample the circumference of the circular cylinder. It will give us a point of the cylinder. - Project it to sphere.(see figure below for Sphere with R = 1)

You get the sample on the sphere. Given below is the algorithm for sampling a point (x,y,x) on a sphere of unit radius, with center at (0,0,0): ~~~ //Choose a height sample uniformly distributed in [-1,1] and set it to z. z = random(&seed0, &seed1)*2-1 //Choose φ uniformly distributed on [0, 2*π). φ = 2*π*random(&seed0, &seed1) r = √(1-z*z) x = r * cosφ y = r * sinφ ~~~ To get a sample on a sphere of radius R and located at (x_c,y_c,z_c) you must use the following additional steps. ~~~ x = x_c + R*x y = y_c + R*y z = z_c + R*z ~~~ `35:T9e5,md`### Sampling for direction on a Unit Hemisphere around vector N Let us assume we know how to sample a direction d (x_d,y_d,z_d) on a unit hemisphere around Z-axis. Then all we have to do is transform d by a rotation matrix such that the rotation aligns Z-axis with the vector N. Instead of looking for rotation we will use our knowledge of coordinate transformation to compute this transformation matrix. That means we will set up a UVW coordinate axes such that W=N, and U,V are orthogonal unit vector perpendicular to N. Then the transformation matrix will be: [U V W]. We will will apply this matrix to d and compute the required direction vector. i.e. ~~~ dsample = transpose([U V W]) * d = [dot(d*U) , dot(d*V), dot(d*W)]. ~~~ Given that we know W, all we need is two orthogonal unit vectors U,V that are perpendicular to W. Here is an algorithm. Assume N be normalized vector. ~~~ tempVector=N if (|x_N| < |y_N| && |x_N| < |z_N|) x_tempVector=1; else if (|y_N) < |x_N| && |y_N| < |z_N|) y_tempVector=1; else z_tempVector=1; U = normalize(cross(tempVector,N)) V = cross(N,U) W = N ~~~ Note that any vector different from vector N can serve as tempVector. You may even do the following: ~~~ if (|x_N| > .1f) tempVector = (0.f, 1.f, 0.f); else tempVector = (1.f, 0.f, 0.f); ~~~ #### Sampling direction uniformly on a Unit Hemisphere Earlier we discussed a method for sampling point on a unit sphere whose center is at (0,0,0). That same method may be used with a minor change. The algorithm is as follows: ~~~ //Choose a height sample uniformly distributed in [0,1] and set it to z_d. z = random(&seed0, &seed1) // Choose φ uniformly distributed on [0, 2*pi). φ = 2*π*random(&seed0, &seed1) r = √(1-z*z) xd = r * cosφ yd = r * sinφ ~~~ ### Sampling direction according to cosine function on a Unit Hemisphere around Z-axis This method is based on the fact that projection of the surface of a unit hemisphere on its base plane is a unit disc. We know what role cosine function plays in projections. Earlier on this page, we described a disc sampling method of sampling hemisphere. So if we uniformly sample that disc and unproject the sample points then we will get points distributed over a hemisphere according to cosine distribution. The algorithm is as follows: ~~~ // Sample a point on the disc. Any of the earlier mentioned method may be used. φ = 2*π* random(&seed0, &seed1) r = √random(&seed0, &seed1) x_d = r*cosφ y_d = r*sinφ //Unproject the point: only Z component changes z_d = √(1-r*r) ~~~ `36:[{"id":"31279eb8ad6436a099009458ea7ff18fdeb182227a74099397d6f0a345af6634bf76854cc2eb93f13b0e7a66c2cea1998bbfc49dba57419fa380a1090205fc63","href":"https://static.observableusercontent.com/files/31279eb8ad6436a099009458ea7ff18fdeb182227a74099397d6f0a345af6634bf76854cc2eb93f13b0e7a66c2cea1998bbfc49dba57419fa380a1090205fc63","name":"cylinderPoint2sphere.png","type":"image/png","lastModified":1614031646409.609,"size":385015}] 2e:{"title":"Sampling","cells":[{"id":0,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":2,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":18,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":28,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":72,"value":"A Simplification of the code from a [graphics blog](http://psgraphics.blogspot.com/2011/01/improved-code-for-concentric-map.html) and [a code from github page](https://github.com/erich666/jgt-code/blob/master/Volume_02/Number_3/Shirley1997/disk.cc).\n","pinned":false,"hidden":false,"mode":"md"},{"id":74,"value":"~~~\nif (sx == 0.0 && sy == 0.0) {// Handle degeneracy at the origin\n x = 0;\n y = 0;\n}\nelse{\n if (sx*sx > sy*sy){\n r = sx;\n θ = (π / 4)*(sy/sx);\n }\n else{\n r = sy;\n θ = (π / 2)-(π / 4)*(sx/sy)\n }\n x = r * cos(θ);\n y = r * sin(θ);\n}\n~~~","pinned":false,"hidden":false,"mode":"md"},{"id":30,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":32,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W36","stdlib":"1","titling":"auto"}