Solid angle of triangle

[tuomlarsen]

Hello,

I'm trying to understand how a path tracer works and in the process I thought I should try to extend smallpt http://kevinbeason.com/smallpt/explicit.cpp to support triangles.

In the explicit light sampling part, there it sums all the unshadowed lights:

Code: Select all

double omega = 2*M_PI*(1-cos_a_max);
e = e + f.mult(s.e*l.dot(nl)*omega)*M_1_PI;

So the task is to compute "omega" for triangles. I implemented the method described at http://en.wikipedia.org/wiki/Solid_angle#Tetrahedron and it seems to work.

But there is one thing which bothers me. Solid angle is area divided by squared distance, right?

Code: Select all

omega = A / r * r

By cosine law, I should be able to "scale" the area by the dot product of light ray direction and the light triangle normal. (Btw., is this called "projected solid angle"?) In the above code, it would read:

Code: Select all

scale = l.dot(normal_of_s)

I tried it, and it seems to work, too.

So my questions is, or rather, my questions are:

- Why does it work? I mean, obviously, the two methods are different and the second one is much more simple - is the second method just an approximation which works only if, say, the solid angles are small enough?

- If it works, why to use the first method at all? I think I'm totally missing something, I would be thankful if someone could point me to the right direction.

Thanks in advance!

[dr_eck]

Thanks for posting this basic question on the "My First" forum. As described in the Wikipedia article, solid angle is defined as the area of the triangle projected onto a unit sphere with its center at the observation point (the point from which your "r" is measured).

As a counterexample to your formula, omega = A/r*r, what if r = 1 and A = 1e6? The correct solid angle in this case (assuming a planar area perpendicular to r) would be approximately pi, not 1e6. Your formula is merely an approximation that works well for A << r*r.

Regarding your question about projected solid angle: Yes, if you are going to use your approximate formula, a triangle seen edge-on will have 0 solid angle, so it is necessary to use the dot product (i.e. the projected area), not just the area. If you want to get really picky, "projected solid angle" is redundant, because solid angle is always based on the projected area.

[tuomlarsen]

Thank you very much for the reply, much appreciated!

One more question, though:

As a counterexample to your formula, omega = A/r*r, what if r = 1 and A = 1e6? The correct solid angle in this case (assuming a planar area perpendicular to r) would be approximately pi, not 1e6.

Solid angle is the area "cut out" from unit sphere surface divided by 4*pi (as the total surface area of unit sphere). So if r = 1 and A = 1e6, the omega should be 1e6 / (4*pi), rather than pi, no?

PS: Now reading Wikipedia, there is the formula: omega = 2 * pi (1 - cos theta). I would assume that since the area is tiny, the angle should also be tiny, because cos(tiny angle) is 1, so omega = 2 * pi (1 - 1)..?

[beason]

I believe I wrote some code to compute the exact projected solid angle of a triangle clipped to the horizon.
This is from http://svn.kevinbeason.com/pane/trunk/PaneDirect.cxx

Code: Select all

      // compute projected solid area using Stoke's thm
      // from Improving Radiosity Solutions
      // through the Use of Analytically Determined
      // Form Factors by Baum, Rushmeier, and Winget
      // (Eq. 9 on pg. 6 (or "330"))

      SbVec3f pts[4];
      int ptscnt = 0;

      // clip triangle against horizon
      static const int nextIdx[] = { 1, 2, 0 };  // j = (i+1)%3
      for (int i=0; i<3; i++){
        int j = nextIdx[i];
        if (dp[i]>=0.0)      // check vertex (don't clip pts on horizon)
          pts[ptscnt++] = vert[i];
        if (dp[i]*dp[j]<0.0) // clip crossing edge against horizon
          pts[ptscnt++] = (vert[i]*(-dp[j]) + vert[j]*dp[i])/(dp[i]-dp[j]);
      }

      // guaranteed to be 3-4 points
      // integrate over edges
      float projarea = 0.0;
      for (int i=0; i<ptscnt; i++){
        int j = (i+1)%ptscnt;
        SbVec3f r1 = pts[i]-r.x;
        SbVec3f r2 = pts[j]-r.x;
        SbVec3f tau = r1.cross(r2);
        r1.normalize();
        r2.normalize();
        float dotp = clamp(r1.dot(r2), -1.f, 1.f);
	
        float gamma = acos(dotp);
        if (!finite(gamma)){
          cerr << "gamma not finite" << endl;
          exit(1);
        }
      
        tau.normalize();
        tau *= gamma;
        projarea -= ng.dot(tau);
	
        if (!finite(projarea)){
          cerr << "problem with projarea\n";
          exit(1);
        }
	
      }

      //projarea /= 2.0*M_PI;  // This would give integral of 1...
      //projarea *= M_PI;      // ...but we want pi
      projarea *= 0.5;         // so I think 1/2 is the normalization we want for the integral
 

[dr_eck]

Thank you very much for the reply, much appreciated!

One more question, though:

As a counterexample to your formula, omega = A/r*r, what if r = 1 and A = 1e6? The correct solid angle in this case (assuming a planar area perpendicular to r) would be approximately pi, not 1e6.

Solid angle is the area "cut out" from unit sphere surface divided by 4*pi (as the total surface area of unit sphere). So if r = 1 and A = 1e6, the omega should be 1e6 / (4*pi), rather than pi, no?

PS: Now reading Wikipedia, there is the formula: omega = 2 * pi (1 - cos theta). I would assume that since the area is tiny, the angle should also be tiny, because cos(tiny angle) is 1, so omega = 2 * pi (1 - 1)..?

Oops! I should have said that the solid angle is 2*pi, not pi. Also, 1e6 = 1000000, which is not tiny. Perhaps you were thinking of
1e-6 (= 0.0000001)?

Edit: (Gotta try the new MathJax support!)
\[
\Omega = 2 \pi \left( 1-\cos \theta \right)
\]
I'm not comfortable with your term "cut out", probably because my many years of school bias me. But if you can figure out a way to cut out 1000000 square meters of cloth from a piece that is only 4*pi square meters, I think you will become very rich.

Here is a better description of my thinking: The surface area of a unit sphere is 4*pi, as you stated. The 1e6 area would project onto nearly half of the sphere, so the solid angle is half of the area of the sphere, or 2*pi.

One way to picture projection is to imagine an infinite sphere concentric with the unit sphere. The infinite sphere radiates light inward toward the unit sphere, but all of the light rays are normal (perpendicular) to the surface of the sphere. The projection of any area onto the unit sphere corresponds exactly to the shadow that would be generated. In the case of the very large area, roughly half of the unit sphere would be shadowed.

Another way of thinking about this projection is from the inside out. If you trace rays at random angles from the center of the sphere, roughly half of them would hit the 1e6 area, so the solid angle of the area must be half that of the unit sphere.

Using the last formula you give, omega = 2 * pi * (1 - cos(pi/2)) = 2*pi*(1-0) = 2*pi

Another MathJax Edit:
\[
\Omega = 2 \pi \left( 1 - \cos \left( \pi / 2 \right) \right) = 2 \pi \left( 1 - 0 \right) = 2 \pi
\]

Ooh, I like MathJax!

[tuomlarsen]

dr_eck,

thank you very much! I now see where the simpler method breaks. And stupid me, 1e6 is not tiny at all...

beason,

thanks for the posted code, I found the whole file interesting, especially the part where you choose between approximate and accurate sampling. Also, I (and I'm sure many others) cannot thank you enough for writing smallpt, it helps a lot to see how the basics fit together.