Monday, June 4, 2018

Sticking a thin lens in a ray tracer

I need to stick a " ideal thin lens" in my ray tracer.   Rather than a real lens with some material, coating, and shape, it's an idealized version like we get in Physics1.

The thin lens has some basic properties that have some redundancy but they are the ones I remember and deploy when needed:
  1. a ray through the lens center is not bent
  2. all rays through a point that then pass through the lens converge at some other point
  3. a ray through the focal point (a point on the optical axis at distance from lens f, the focal length of the lens) will be parallel to the optical axis after being bent by the lens
  4. all rays from a point a distance A from the lens will converge at a distance B on the other side of the lens and obey the thin lens law: 1/A + 1/B = 1/f
 Here's a sketch of those rules:

So if I have a (purple) ray with origin a that hits the lens at point m, how does it bend?   It bends towards point b no matter what the ray direction is.   So the new ray is:

p(t) = m + t (b-m)

 So what is point b?

It is in the direction of point c, but extended by some distance.  We know the center of the lens c, so we can use the ratios of the segments to extend to that:

b = a + (c-a)  (B+A) / A?

So what is (B+A) /A? 

We know 1/A + 1/B = 1/f, so B = 1/(1/f-1/A).   So the point b is:

b = a + (c-a) (1/(1/f-1/A) + A)/A =
b = a + (c-a) (1/(A/f - 1) + 1) =
b = a + (c-a) (A/(A - f))

OK lets try a couple of trivial cases.   What if  A = 2f?  

b = a + (c-a) (2f/(2f-f)) = b = a + 2(c-a) 

That looks right (symmetric case-- A = B there) 

So final answer, given a ray with origin a that hits a lens with center c and focal length f at point m, the refracted ray is:

p(t) = m + t( a + (c-a) (A/(A - f)) - m)

There is a catch.   What if B < 0?   This happens when A < f.  Address that case when it comes up :)

Thursday, May 31, 2018

Generating uniform Random rays that hit an axis aligned box

For some tests you want the set of "all" rays that hit a box.   If you want to stratify this is somewhat involved (and I don't know that I have done it nor seen it).   Chris Wyman and I did a jcgt paper on doing this stratified in a square in 2D.   But often stratification isn't needed.   When in doubt I never do it-- the software is always easier to deal with un-stratified, and as soon as dimension gets high most people dont bother because of the Curse of Dimensionality.

We would like to do this with as little math as possible.   First lets consider any side of the box.   (this would apply to any convex polyhedron if we wanted something more general).    If the box in embedded in all possible uniform rays, any ray that hits the box will enter at exactly one point on the surface of the box, and all points are equally likely.   So our first job is to pick a uniform point on the surface.   We can use the technique Greg Turk used to seed points for texture generation on a triangular mesh:

Probability of each side with box of side lengths X, Y, Z is side area of total area.   The side areas are:

XY, YZ, ZX (2 of each)

We can do cumulative area and stuff it in an array of length 6:

c_area[0] = XY
c_area[1] = area[0] + XY
c_area[2] = area[1] + YZ
c_area[3] = area[2] + YZ
c_area[4] = area[3] + ZX
c_area[5] = area[4] + ZX 

Now normalize it so it is a cumulative fraction:

for (int i = 0; i < 6; i++) 
    c_area[i] /= c_area[5]

No take a uniform random real r() in [0,1)

int candidate = 0;
float ra = r();
while (c_area[candidate] < ra) candidate++;

Now candidate is the index to the side.

Let's say the side is in the xy plane and x goes from 0 to X and y goes from 0 to Y.   Now pick a uniform random point on that side:

vec3 ray_entry(X*r(), Y*r(), 0);

Now we have a ray origin.   What is the ray direction?  It is not uniform in all directions.   These are the rays that hit the side.   So the density is proportional to the cosine to the normal-- so they are Lambertian!   This is not obvious.   I will punt on justifying that for now.

So for the xy plane one above, the normal is +z and the ray direction is a uniform random point on a disk projected onto the hemisphere:

float radius = sqrt(r());
float theta = 2*M_PI*r();
float x = radius*cos(theta);
float y = radius*sin(theta);
float z = sqrt(1-x*x-y*y);
ray_direction = vec3(x,y,z);

Now we need that for each of the six sides.

We could probably find symmetries to have 3 cases, or maybe even a loop, but I personally would probably not bother because me trying to be clever usually ends poorly...

Wednesday, March 14, 2018

Egyptian estimates of PI

I saw a neat tweet on how the estimate the Egyptians used for PI.

This is all my speculation, and maybe a math history buff can enlighten me, but the D^2 dependence they should discover pretty naturally.   Including the constant before squaring is, I would argue, just as natural as having it outside the parentheses, so let's go with that for now.   So was there a nearby better fraction?   How well did the Egyptians do?   A brute force program should tell us.

We will use the ancient programming language C:

int main() {
    double min_error = 10;
    for (int denom = 1; denom < 10000; denom++) {
        for (int num = 1; num < denom; num++) {
            double approx = 2*double(num)/double(denom);
            approx = approx*approx;
            double error2 = M_PI-approx;
            error2 = error2*error2;
            if (error2 < min_error) {
                min_error = error2;
                printf("%d/%d %f\n", num, denom, 4*float(num*num)/float(denom*denom));

This produces output:

1/2 1.000000
2/3 1.777778
3/4 2.250000
4/5 2.560000
5/6 2.777778
6/7 2.938776
7/8 3.062500
8/9 3.160494
23/26 3.130177
31/35 3.137959
39/44 3.142562
109/123 3.141252
148/167 3.141597
4401/4966 3.141588
4549/5133 3.141589
4697/5300 3.141589
4845/5467 3.141589
4993/5634 3.141589
5141/5801 3.141590
5289/5968 3.141590
5437/6135 3.141590
5585/6302 3.141590
5733/6469 3.141590
5881/6636 3.141591
6029/6803 3.141591
6177/6970 3.141591
6325/7137 3.141591
6473/7304 3.141591
6621/7471 3.141591
6769/7638 3.141591
6917/7805 3.141592
7065/7972 3.141592
7213/8139 3.141592
7361/8306 3.141592
7509/8473 3.141592
7657/8640 3.141592
7805/8807 3.141592
7953/8974 3.141593
8101/9141 3.141593
8249/9308 3.141593
8397/9475 3.141593
8545/9642 3.141593

So 7/8 was already pretty good, and you need to get to 23/26 before you do any better!   I'd say the Egyptians did extremely well.

What if they had put the constants outside the parens?   How well could they have done?   We can change two of the lines above to:

double approx = 4*double(num)/double(denom);//approx = approx*approx;

and the printf to:

printf("%d/%d %f\n", num, denom, 4*float(num)/float(denom));

And we get:

1/2 2.000000
2/3 2.666667
3/4 3.000000
4/5 3.200000
7/9 3.111111
11/14 3.142857
95/121 3.140496
106/135 3.140741
117/149 3.140940
128/163 3.141104
139/177 3.141243
150/191 3.141361
161/205 3.141464
172/219 3.141552
183/233 3.141631
355/452 3.141593

So 7/9 is not bad!  And 11/14 even better.   So no clear winner here on whether the rational constant should be inside the parens or not.