Added projection ray blaster & options for periodic blaster/scene

Author: langmm

hothouse/blaster.py

@@ -32,6 +32,34 @@ class RayBlaster(traitlets.HasTraits):
     intensity = traitlets.CFloat(1.0)
     diffuse_intensity = traitlets.CFloat(0.0)
     multibounce = traitlets.CBool(False)
+    period = traittypes.Array(
+        np.zeros((3,), "f4")
+    ).valid(check_dtype("f4"), check_shape(3))
+    periodic_direction = traittypes.Array(np.array([
+        [0.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]
+    ], "f4")).valid(check_dtype("f4"), check_shape(3, 3))
+    periodic_count = traittypes.Array(np.array([1, 1, 1], "i4")).valid(
+        check_dtype("i4"), check_shape(3))
+
+    @classmethod
+    def from_scene(cls, scene, **kwargs):
+        r"""Create a blaster for a scene by setting missing parameters
+        that optimize the coverage of the scene.
+
+        Args:
+            scene (hothouse.scene.Scene): Scene to customize the blaster
+                for.
+
+        Returns:
+            RayBlaster: Created blaster.
+
+        """
+        raise NotImplementedError
+
+    @cached_property
+    def is_periodic(self):
+        r"""bool: True if the blaster will be replicated periodically."""
+        return np.any(self.period > 0)
 
     @property
     def ray_intensity(self):
@@ -43,20 +71,38 @@ def cast_once(self, scene, verbose_output=False,
         if self.multibounce:
             callback_handler = RayCollisionMultiBounce(
                 self.origins.shape[0], 10,
-                [c.transmittance for c in scene.components],
-                [c.reflectance for c in scene.components])
+                scene.transmittance, scene.reflectance)
         else:
             callback_handler = None
+        dists = None
+        if self.is_periodic:
+            dists = np.empty(self.origins.shape[0], 'float32')
+            dists.fill(1e37)
         output = scene.embree_scene.run(
             self.origins,
             self.directions,
+            dists=dists,
             query=query_type._value_,
             output=verbose_output,
             callback_handler=callback_handler
         )
+        if self.is_periodic:
+            from hothouse.scene import PeriodicScene
+            shifts = PeriodicScene.get_periodic_shifts(
+                self.period, self.periodic_direction,
+                self.periodic_count)
+            for shift in shifts:
+                output = scene.embree_scene.run(
+                    self.origins + shift,
+                    self.directions,
+                    dists=dists,
+                    query=query_type._value_,
+                    output=verbose_output,
+                    callback_handler=callback_handler
+                )
         if self.multibounce and isinstance(output, dict):
             output['bounces'] = callback_handler.bounces
-        return output
+        return scene.post_cast(output)
 
     def compute_distance(self, scene):
         output = self.cast_once(
@@ -159,6 +205,13 @@ class SunRayBlaster(OrthographicRayBlaster):
     scene_limits = traittypes.Array(None, allow_none=True).valid(
         check_shape(None, 3), check_dtype("f4"))
 
+    @classmethod
+    def get_solar_direction(cls, latitude, longitude, date, up, north):
+        instance = cls(latitude=latitude, longitude=longitude,
+                       date=date, zenith=(10 * up),
+                       ground=np.zeros((3,), "f4"), north=north)
+        return instance.forward
+
     @traitlets.default("_solpos_info")
     def _solpos_info_default(self):
         return pvlib.solarposition.get_solarposition(
@@ -286,7 +339,56 @@ def _up_default(self):
 
 
 class ProjectionRayBlaster(RayBlaster):
-    pass
+    fov_width = traitlets.CFloat(90.0)
+    fov_height = traitlets.CFloat(90.0)
+    center = traittypes.Array().valid(check_dtype("f4"), check_shape(3))
+    forward = traittypes.Array().valid(check_dtype("f4"), check_shape(3))
+    up = traittypes.Array().valid(check_dtype("f4"), check_shape(3))
+    east = traittypes.Array().valid(check_dtype("f4"), check_shape(3))
+    width = traitlets.CFloat(1.0)
+    height = traitlets.CFloat(1.0)
+    nx = traitlets.CInt(512)
+    ny = traitlets.CInt(512)
+
+    @cached_property
+    def camera_distance(self):
+        r"""float: Distance of the camera from the image plane."""
+        return (
+            (self.width / 2.0)
+            / np.tan(np.radians(self.fov_width / 2.0))
+        )
+
+    @cached_property
+    def camera_origin(self):
+        r"""np.ndarray: Position of the camera."""
+        return self.center - self.camera_distance * self.forward
+
+    @traitlets.default("east")
+    def _default_east(self):
+        return np.cross(self.forward, self.up)
+
+    @traitlets.default("directions")
+    def _default_directions(self):
+        self._directions = self.origins - self.camera_origin
+        norms = np.linalg.norm(self.origins, axis=1)
+        for i in range(3):
+            self._directions[:, i] /= norms
+        return self._directions.view()
+
+    @traitlets.default("origins")
+    def _default_origins(self):
+        # here origin is not the center, but the bottom left
+        self._origins = np.zeros((self.nx, self.ny, 3), dtype="f4")
+        offset_x, offset_y = np.mgrid[
+            (-self.width / 2):(self.width / 2):(self.nx * 1j),
+            (-self.height / 2):(self.height / 2):(self.ny * 1j),
+        ]
+        self._origins[:] = (
+            self.center
+            + offset_x[..., None] * self.east
+            + offset_y[..., None] * self.up
+        )
+        return self._origins.view().reshape((self.nx * self.ny, 3))
 
 
 class SphericalRayBlaster(RayBlaster):

hothouse/scene.py

@@ -2,6 +2,7 @@
 import traitlets
 import pythreejs
 import numpy as np
+import functools
 import pvlib
 # from IPython.core.display import display
 
@@ -13,10 +14,13 @@
 from pyembree import rtcore_scene as rtcs
 from pyembree.mesh_construction import TriangleMesh
 
+cached_property = getattr(functools, "cached_property", property)
+
 
 class Scene(traitlets.HasTraits):
-    
-    ground = traittypes.Array(np.array([0.0, 0.0, 0.0], "f4")).valid(check_dtype("f4"), check_shape(3))
+
+    ground = traittypes.Array(np.array([0.0, 0.0, 0.0], "f4")).valid(
+        check_dtype("f4"), check_shape(3))
     up = traittypes.Array(np.array([0.0, 0.0, 1.0], "f4")).valid(
         check_dtype("f4"), check_shape(3))
     north = traittypes.Array(np.array([0.0, 1.0, 0.0], "f4")).valid(
@@ -29,9 +33,14 @@ class Scene(traitlets.HasTraits):
     # TODO: Add surface for ground so that reflection from ground
     # is taken into account
 
+    def post_cast(self, output):
+        return output
+
     def add_component(self, component):
-        self.components = self.components + [component]  # Force traitlet update
-        self.meshes.append(TriangleMesh(self.embree_scene, component.triangles))
+        # Force traitlet update
+        self.components = self.components + [component]
+        self.meshes.append(
+            TriangleMesh(self.embree_scene, component.triangles))
 
     def compute_hit_count(self, blaster):
         output = blaster.compute_count(self)
@@ -43,6 +52,30 @@ def compute_hit_count(self, blaster):
             )
         return component_counts
 
+    @property
+    def transmittance(self):
+        return [c.transmittance for c in self.components]
+
+    @property
+    def reflectance(self):
+        return [c.reflectance for c in self.components]
+
+    @cached_property
+    def limits(self):
+        r"""np.ndarray: Positions of corners of a box containing all
+        points in the scene."""
+        mins = []
+        maxs = []
+        for c in self.components:
+            mins.append(np.min(c.vertices, axis=0))
+            maxs.append(np.max(c.vertices, axis=0))
+        mins = np.min(np.vstack(mins), axis=0)
+        maxs = np.max(np.vstack(maxs), axis=0)
+        limits = np.vstack([mins, maxs])
+        xx, yy, zz = np.meshgrid(limits[:, 0], limits[:, 1], limits[:, 2])
+        limits = np.vstack([xx.flatten(), yy.flatten(), zz.flatten()]).T
+        return limits
+
     def get_sun_blaster(self, latitude, longitude, date,
                         direct_ppfd=1.0, diffuse_ppfd=1.0, **kwargs):
         r"""Get a sun blaster that is adjusted for this scene so that
@@ -72,23 +105,14 @@ def get_sun_blaster(self, latitude, longitude, date,
         # TODO: Calculate direct/diffuse ppfd from lat/long/date
         # using pvi if not provided
         max_distance2 = 0.0
-        mins = []
-        maxs = []
         for c in self.components:
-            mins.append(np.min(c.vertices, axis=0))
-            maxs.append(np.max(c.vertices, axis=0))
             max_distance2 = max(
                 max_distance2,
                 np.max(np.sum((c.vertices-self.ground)**2, axis=1)))
-        mins = np.min(np.vstack(mins), axis=0)
-        maxs = np.max(np.vstack(maxs), axis=0)
-        limits = np.vstack([mins, maxs])
-        xx, yy, zz = np.meshgrid(limits[:, 0], limits[:, 1], limits[:, 2])
-        limits = np.vstack([xx.flatten(), yy.flatten(), zz.flatten()]).T
         max_distance = np.sqrt(max_distance2)
         kwargs.setdefault('zenith', self.up * max_distance + self.ground)
         kwargs.setdefault('diffuse_intensity', diffuse_ppfd)
-        kwargs.setdefault('scene_limits', limits)
+        kwargs.setdefault('scene_limits', self.limits)
         blaster = SunRayBlaster(latitude=latitude,
                                 longitude=longitude, date=date,
                                 ground=self.ground, north=self.north,
@@ -143,7 +167,7 @@ def update(frame):
             o[o <= 0] = np.nan
             img.set_data(o.reshape((camera.ny, camera.nx), order='F'))
             return img,
-            
+
         dates = pd.date_range(t_start, t_stop, periods=n_step)
         ani = FuncAnimation(fig, update, frames=list(dates), blit=False)
         if fname is None:
@@ -226,14 +250,15 @@ def _calc_incident_power(self, ray_dir, norm, area, any_direction=True):
             if any_direction:
                 aoi[aoi > np.pi/2] -= np.pi
             else:
-                aoi[aoi > np.pi/2] = np.pi  # No contribution
+                aoi[aoi > np.pi/2] = np.pi / 2  # No contribution
         else:
             if aoi > np.pi/2:
                 if any_direction:
                     aoi -= np.pi
                 else:
-                    aoi = np.pi
-        return np.cos(aoi) / area
+                    aoi = np.pi / 2
+        out = np.cos(aoi) / area
+        return out
 
     def _accumulate_hits(self, component_fd, primID, geomID,
                          ray_dir, ray_intensity, diffuse_intensity,
@@ -270,13 +295,17 @@ def _accumulate_hits(self, component_fd, primID, geomID,
                 tilt = np.arccos(
                     np.dot(norms, self.up)
                     / (2.0 * areas * np.linalg.norm(self.up)))
-                component_fd[ci] += pvlib.irradiance.isotropic(
+                component_diffuse = pvlib.irradiance.isotropic(
                     np.degrees(tilt), diffuse_intensity)
+                component_fd[ci] += component_diffuse
+            # assert not any(component_fd[ci] == 0)
 
     def _ipython_display_(self):
-        # This needs to actually display, which is not the same as returning a display.
+        # This needs to actually display, which is not the same as
+        # returning a display.
         cam = pythreejs.PerspectiveCamera(
-            position=[25, 35, 100], fov=20, children=[pythreejs.AmbientLight()],
+            position=[25, 35, 100], fov=20,
+            children=[pythreejs.AmbientLight()],
         )
         children = [cam, pythreejs.AmbientLight(color="#dddddd")]
         material = pythreejs.MeshBasicMaterial(
@@ -301,3 +330,102 @@ def _ipython_display_(self):
         )
 
         return rendererCube
+
+
+class PeriodicScene(Scene):
+
+    period = traittypes.Array(
+        np.zeros((3,), "f4")
+    ).valid(check_dtype("f4"), check_shape(3))
+    direction = traittypes.Array(np.array([
+        [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]
+    ], "f4")).valid(check_dtype("f4"), check_shape(3, 3))
+    count = traittypes.Array(np.array([1, 1, 1], "i4")).valid(
+        check_dtype("i4"), check_shape(3))
+    buffer_as_primary = traitlets.Bool(False)
+
+    def __init__(self, *args, **kwargs):
+        self._buffer_meshes = []
+        self._buffer_component_map = {}
+        self._buffer_component_triangles = []
+        self._embree_scene = kwargs.get('embree_scene',
+                                        rtcs.EmbreeScene())
+        super(PeriodicScene, self).__init__(*args, **kwargs)
+
+    @property
+    def transmittance(self):
+        out = super(PeriodicScene, self).transmittance
+        self.add_component_buffers()
+        for i, v in self._buffer_component_map.items():
+            out += [self.components[i].transmittance for _ in v]
+        return out
+
+    @property
+    def reflectance(self):
+        out = super(PeriodicScene, self).reflectance
+        self.add_component_buffers()
+        for i, v in self._buffer_component_map.items():
+            out += [self.components[i].reflectance for _ in v]
+        return out
+
+    @classmethod
+    def get_periodic_shifts(cls, period, direction, count):
+        import itertools
+        shifts = []
+        opts = []
+        for axis in range(3):
+            if period[axis] == 0:
+                opts.append([0])
+                continue
+            else:
+                opts.append(list(range(-count[axis], count[axis] + 1)))
+        for xbuffer, ybuffer, zbuffer in itertools.product(*opts):
+            if xbuffer == 0 and ybuffer == 0 and zbuffer == 0:
+                continue
+
+            def _shift(ibuffer, axis):
+                return ibuffer * period[axis] * direction[axis, :]
+
+            shifts.append(_shift(xbuffer, 0)
+                          + _shift(ybuffer, 1)
+                          + _shift(zbuffer, 2))
+        return np.vstack(shifts)
+
+    def post_cast(self, output):
+        if self.buffer_as_primary:
+            for ci, v in self._buffer_component_map.items():
+                for ci_per in v:
+                    output["geomID"][output["geomID"] == ci_per] = ci
+            for irange, jrange in self._buffer_component_triangles[::-1]:
+                output["primID"][output["primID"] >= jrange.start] -= (
+                    jrange.start - irange.start)
+        return output
+
+    def add_component_buffers(self):
+        if self._buffer_meshes:
+            return
+        shifts = self.get_periodic_shifts(self.period, self.direction,
+                                          self.count)
+        j = len(self.components)
+        icount = 0
+        jcount = sum([component.triangles.shape[0]
+                      for component in self.components])
+        for i, component in enumerate(self.components):
+            irange = range(icount, icount + component.triangles.shape[0])
+            self._buffer_component_map[i] = []
+            for shift in shifts:
+                jrange = range(jcount, jcount + component.triangles.shape[0])
+                self._buffer_meshes.append(
+                    TriangleMesh(self._embree_scene,
+                                 component.triangles + shift)
+                )
+                self._buffer_component_map[i].append(j)
+                self._buffer_component_triangles.append((irange, jrange))
+                j += 1
+                jcount += component.triangles.shape[0]
+            icount += component.triangles.shape[0]
+
+    @property
+    def embree_scene(self):
+        self.add_component_buffers()
+        return self._embree_scene