broken D.ndim == 2 option in shape.py

[steo85it]

I think that while updating shape.py with the is_occluded function, the option for an extended light source (D.ndim == 2, with D the direction between the vector array to the "Sun") got broken.

The current _is_occluded doesn't work in either cgal or embree version. For the cgal version, the issue is visible at

https://github.com/sampotter/python-flux/blob/700965ad79a8f387869459da433f63ad2629fe32/src/flux/shape.py#L263-L267

where D[i] cycles on the number of shape-model's facets instead of the number of "Sun facets/vectors" (in my case, D.shape=(100,3)).

For reference, see the old function get_direct_irradiance in the attached file (for a working embree version).
shape.txt

Any help is welcome @sampotter

[sampotter]

Could you send me a short script which I can use to make this bug occur?

[steo85it]

it took a moment (the original mesh set-up in Haworth.py is suuuuuuper slow...), but I should have set up a "haworth-like" example using an extended Sun (some stuff might be off, physically, but at least I get the same error message) - please see related PR

[sampotter]

Got it, I'm taking a look now...

[sampotter]

The way _is_occluded is implemented now is wrong for D.ndim == 2. In this case, is_occluded should return a len(I) x D.shape[0] boolean matrix. Try this:

def _is_occluded(self, I, D):
    m = len(I)
    if D.ndim == 1:
        occluded = np.empty(m, dtype=np.bool_)
        for p, i in enumerate(I):
            occluded[p] = self.aabb.ray_from_centroid_is_occluded(i, D)
    elif D.ndim == 2:
        n = D.shape[0]
        assert D.shape[1] == 3
        occluded = np.empty((m, n), dtype=np.bool_)
        for (p, i), (q, d) in zip(enumerate(I), enumerate(D)):
            occluded[p, q] = self.aabb.ray_from_centroid_is_occluded(i, d)
    return occluded

I ran your example with this code but I hit another error:

Traceback (most recent call last):
  File "/Users/sfp/Dropbox/Research/ThermalModeling/python-flux-steo85it-extsun/examples/lunar_south_pole/haworth_extsun.py", line 113, in <module>
    E = shape_model.get_direct_irradiance(F0 / 100., sundirs)
  File "/opt/local/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/python_flux-0.0.0-py3.9-macosx-12.0-x86_64.egg/flux/shape.py", line 220, in get_direct_irradiance
    raise ValueError('need Dsun.shape[0] == num_faces')
ValueError: need Dsun.shape[0] == num_faces

So I think there is a slight issue here, because there are two reasonable ways of interpreting the Dsun argument in get_direct_irradiance when Dsun.ndim == 2:

1. Dsun is a collection of sun directions such that Dsun.shape[0] == num_faces, and we assume 1 sun direction per triangle, each of which may be different (e.g. to implement a point sun)
2. Dsun is a collection of multiple sun directions, which should each be used for all triangles (i.e., we want to compute multiple righthand sides for several different "parallel" sun directions)

IMO, in this kind of situation, the best thing to do is just pick a particular semantics for this function and add more functions to handle other types of situations. I don't have a strong opinion about how this is done in this case.

[steo85it]

Thanks, I'm testing it. In the meantime, to clarify, what I'd like to have is the "interpretation n.2", which was implemented in the "old" shape.py that I attached to the issue. We could indeed set up a different function to keep things clearer.

[sampotter]

Please just make a separate call to get_direct_irradiance for each sun direction vector for the time being.

[steo85it]

Ok, no problem, that's also fine although slower. :) I'll keep the issue open, though.
Just for reference, what worked with embree was

n = self.num_faces
...
elif Dsun.ndim == 2:
    m = Dsun.size//3
    ray = embree.Ray1M(m*n)
    for i in range(m):
        ray.org[i*n:(i + 1)*n, :] = self.P + eps[:,np.newaxis]*self.N
    for i, d in enumerate(Dsun):
        ray.dir[i*n:(i + 1)*n, :] = d
    ray.tnear[:] = 0
    ray.tfar[:] = np.inf
    ray.flags[:] = 0
...

which should in principle also work for CGAL (just a matter of managing indexes in the correct way).

[sampotter]

OK, sounds good.

After we submit the paper we can come back to this. I think the right way to address this is to have separate functions for each of the four cases (single sun direction, multiple constant sun directions per triangle, a different single sun direction for each triangle, different multiple sun directions per triangle). They should also have an intuitive and consistent naming convention.

[nschorgh]

I don't necessarily need the Dsun.shape[0] == num_faces option for now, because I decided to work in (x,y,z), so only have one sun direction for the whole globe. It should be easy to restore the original Dsun.ndim == 2 option instead or in addition to that.