5
$\begingroup$

tl;dr I need help finding the code or OSL script for the volume scatter node. Not asking how it works, just a specific, noob-friendly instruction how to find it, or a link to it. Blender is big, and this is my first time, so...

The volume scatter node has an Anisotropy input to control the relative strength of scattering in different directions. When changed from +1 to -1 the scattering changes from more forward to more backward scattering.

It turns out there is a published excerpt from the Blender Cycles Encyclopedia that shows this - I've included two screenshots below.

The (unnormalized) distribution could be as simple as:

I(theta) ~ 1 + A * cos(theta)

where A is the Anisotropy factor, or it could be completely different.

I could build a science experiment in blender to "measure" the intensity distribution as a function of angle. I don't know if there is a "light meter" node yet, but I could do it with cameras, saving rendered images, and analyzing them with NumPy.

But I am wondering if I can just read the math somewhere? Are nodes like this written in OSL, or in C? (I learned about OSL in this helpful answer.) How would I go about finding the code for this node to find how the shape of the distribution is calculated?

enter image description here

enter image description here

| improve this question | |
$\endgroup$
  • 1
    $\begingroup$ Well Blender is open source, so one of the great advantages is you can browse and explore its underlying workings at will. I am no developer and don't know where that node code resides exaclty, but you can browse the Cycles Project Repository yourself or ask a developer where to find these specific features. CHech the "Browse Code" links towards the end. $\endgroup$ – Duarte Farrajota Ramos♦ Apr 4 '17 at 22:50
  • $\begingroup$ @DuarteFarrajotaRamos Thanks! Blender is not small, and I'm not an experienced navigator within large projects. So I'm hoping that with my question I have in fact "... ask(ed) a developer where to find these specific features." with the sentence "How would I go about finding the code for this node..." $\endgroup$ – uhoh Apr 4 '17 at 22:55
  • $\begingroup$ Just stating that because active developers are always busy, and don't often hang about these parts, let alone answer questions. So chances are small that they will naturally see this and you get a proper answer. Maybe try an official mailing list like BF Committers, you may have more luck there. $\endgroup$ – Duarte Farrajota Ramos♦ Apr 4 '17 at 23:03
  • 1
    $\begingroup$ @DuarteFarrajotaRamos OK thanks, This particular question is pretty simple. I think there are a few experienced people left in Blender SE, so It is possible someone here knows, let's see. $\endgroup$ – uhoh Apr 4 '17 at 23:13
  • 1
    $\begingroup$ I'm no developer, but as far as I can tell from looking at it, there is no math in the short bit of code in the osl version of the shader. It can be found here. For the standard shader, it's about the same, at a different location in the source. $\endgroup$ – Timaroberts Apr 7 '17 at 5:07
3
$\begingroup$

@DuarteFarrajotaRamos 's comment turns out to have been a good suggestion - ask someone. I found the IRC channel on one of the linked pages, asked, and got the following link straight away. Exactly what I needed!

https://developer.blender.org/diffusion/B/browse/master/intern/cycles/kernel/closure/volume.h

note: it's the .h file that has the goodies I'm after, not the .c. Line 35 has the math (several lines shown here):

/* HENYEY-GREENSTEIN CLOSURE */

/* Given cosine between rays, return probability density that a photon bounces
 * to that direction. The g parameter controls how different it is from the
 * uniform sphere. g=0 uniform diffuse-like, g=1 close to sharp single ray. */
ccl_device float single_peaked_henyey_greenstein(float cos_theta, float g)
{
  return ((1.0f - g * g) / safe_powf(1.0f + g * g - 2.0f * g * cos_theta, 1.5f)) * (M_1_PI_F * 0.25f);
};

Below are some plots of the shape as a function of g, along with the python script to make the plots them. The Henyey-Greenstein dates back to 1941 in a paper analyzing astronomical scattering of dust. The distribution is not profoundly magic, but it has some roots in science and works nicely, with one parameter g from -1 through 0 to +1 reproducing backwards scattering to isotropic to forward scattering.

If you are interested, here are a number of links with a wide variety of depth and focus. Plots are below that.

http://www.oceanopticsbook.info/view/scattering/the_henyeygreenstein_phase_function

http://adsabs.harvard.edu/abs/1941ApJ....93...70H

http://articles.adsabs.harvard.edu/cgi-bin/nph-iarticle_query?1941ApJ....93...70H&data_type=PDF_HIGH&whole_paper=YES&type=PRINTER&filetype=.pdf

https://www.astro.umd.edu/~jph/HG_note.pdf

http://omlc.org/education/ece532/class3/hg.html

https://en.wikipedia.org/wiki/Monte_Carlo_method_for_photon_transport

http://articles.adsabs.harvard.edu/cgi-bin/nph-iarticle_query?1985A%26A...146...67H&defaultprint=YES&filetype=.pdf

https://www.cs.dartmouth.edu/~wjarosz/publications/dissertation/chapter4.pdf

http://web.cs.wpi.edu/~emmanuel/courses/cs563/S07/talks/Paulo_volumeScattering_wk10_p1.pdf

enter image description here

The function is already available in OSL. This should not be a surprise since it is so fundamental to rendering and scattering simulations of all varieties! This is from v1.7 of the Open Shading Language 1.7 Language Specification - which is outdated. A new and/or better link to 1.8 or whatever is newest is appreciated!

enter image description here

def pu(cos_theta, g):

    top = oneoverfourpi * (1. - g**2)
    bot = (1 + g**2 - 2*g*cos_theta)**1.5

    return top / bot


import numpy as np
import matplotlib.pyplot as plt

halfpi, pi, twopi, fourpi = [f*np.pi for f in [0.5, 1.0, 2.0, 4.0]]
oneoverfourpi = 1. / fourpi

degs, rads = 180/pi, pi/180

theta      = np.linspace(0, pi, 200)
cos_theta, sin_theta = np.cos(theta), np.sin(theta)

dtheta     = theta[1] - theta[0]

g          = np.arange(0.9, -1, -0.15)[:, None]

Pu         = pu(cos_theta, g)

plt.figure()
plt.subplot(3, 1, 1)
for thing in Pu:
    plt.plot(degs*theta, thing)

plt.subplot(3, 1, 2)
for gval, thing in zip(g, Pu):

    if np.abs(gval) < 1E-05:
        plt.plot(degs*theta, thing, '--k')
    else:
        plt.plot(degs*theta, thing)
plt.yscale('log')

plt.subplot(3, 1, 3)
for thing in Pu:
    plt.plot(degs*theta, twopi * sin_theta * dtheta * thing)
plt.show()
| improve this answer | |
$\endgroup$
  • $\begingroup$ @DuarteFarrajotaRamos thank you for your recommendation - I asked! $\endgroup$ – uhoh Apr 7 '17 at 19:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.