[#973] Advancing on _generate_rays()

tofu/data/_class08_generate_rays.py

@@ -11,6 +11,9 @@
 import datastock as ds
 
 
+from ..geom import _GG
+
+
 # ###############################################################
 # ###############################################################
 #                      Main
@@ -212,14 +215,14 @@ def _random(
     min1, max1 = np.min(out_arr[:, 1]), np.max(out_arr[:, 1])
 
     # areas to get safety factor
-    area = doptics['pix_area']
+    area = coll.dobj['camera'][kcam]['dgeom']['pix_area']
     area_max = (max0 - min0) * (max1 - min1)
 
     # start point coords in pixel's plane
     start0, start1 = _seed_pix(
         path=path,
         nrays=nrays,
-        nrays_factor=np.ceil(area_max/area) + 2,
+        nrays_factor=int(np.ceil(area_max/area) + 2),
         min0=min0,
         max0=max0,
         min1=min1,
@@ -228,30 +231,84 @@ def _random(
 
     # shared apertures ?
     pinhole = doptics['pinhole']
+    parallel = coll.dobj['camera'][kcam]['dgeom']['parallel']
 
     # shape camera
     shape_cam = coll.dobj['camera'][kcam]['dgeom']['shape']
 
+    # camera vector
+    # get camera vectors
+    dvect = coll.get_camera_unit_vectors(kcam)
+    lc = ['x', 'y', 'z']
+    if parallel is True:
+        e0i_x, e0i_y, e0i_z = [dvect[f"e0_{kk}"] for kk in lc]
+        e1i_x, e1i_y, e1i_z = [dvect[f"e1_{kk}"] for kk in lc]
+
     # -----------------------------------
     # pinhole camera (shared apertures)
     # -----------------------------------
 
     if pinhole is True:
 
+        # get pixel centers
         cx, cy, cz = coll.get_camera_cents_xyz(kcam)
         cent_cam = np.r_[np.mean(cx), np.mean(cy), np.mean(cz)]
 
-        end0, end1 = _seed_optics(
+        # get end points
+        end0, end1, op_cent, op_e0, op_e1 = _seed_optics(
             coll=coll,
             cent_cam=cent_cam,
             nrays=nrays,
             optics=doptics['optics'],
-            cent_cam=coll.get_camera_cents_xyz(kcam),
         )
 
+        # full versions
+        end0f = np.tile(end0, nrays)
+        end1f = np.tile(end1, nrays)
+
+        # full versions
+        start0f = np.repeat(start0, end0.size)
+        start1f = np.repeat(start1, end0.size)
+
+        # vignetting polygons
+        vignett_poly = [
+            np.array(coll.get_optics_poly(k0, closed=True, add_points=False))
+            for k0 in doptics['optics']
+        ]
+        lnvert = np.array([vv.shape[1] for vv in vignett_poly], dtype=np.int64)
+
+        # loop on pixels
+        ray_orig = np.full((3, nrays*end0.size), np.nan)
+        ray_vdir = np.full((3, nrays*end0.size), np.nan)
         for ind in np.ndindex(shape_cam):
-            pass
 
+            # unit vectors
+            if parallel is not True:
+                e0i_x, e0i_y, e0i_z = [dvect[f"e0_{kk}"][ind] for kk in lc]
+                e1i_x, e1i_y, e1i_z = [dvect[f"e1_{kk}"][ind] for kk in lc]
+
+            # ray_orig
+            ray_orig[0, :] = cx[ind] + start0f * e0i_x + start1f * e1i_x
+            ray_orig[1, :] = cy[ind] + start0f * e0i_y + start1f * e1i_y
+            ray_orig[2, :] = cz[ind] + start0f * e0i_z + start1f * e1i_z
+
+            # ray_vdir, normalized
+            ray_vdir[0, :] = op_cent[0] + end0f * op_e0[0] + end1f * op_e1[0]
+            ray_vdir[1, :] = op_cent[1] + end0f * op_e0[1] + end1f * op_e1[1]
+            ray_vdir[2, :] = op_cent[2] + end0f * op_e0[2] + end1f * op_e1[2]
+            ray_norm = np.sqrt(np.sum(ray_vdir**2, axis=0))
+            ray_vdir[:] = ray_vdir / ray_norm
+
+            # vignetting (npoly, nrays)
+            iok = _GG.vignetting(
+                ray_orig,
+                ray_vdir,
+                vignett_poly,
+                lnvert,
+                num_threads=16,
+            )
+
+            print(ind, iok.sum(), iok.size)
 
 
     # ---------------
@@ -332,15 +389,15 @@ def _seed_optics(
     # seed
     # -----------------------------
 
-    out0, out1 = coll.get_optics_outline(kop)
+    out0, out1 = coll.get_optics_outline(kop, add_points=False, closed=False)
 
     path = mpath.Path(np.array([out0, out1]).T)
 
     min0, max0 = np.min(out0), np.max(out0)
     min1, max1 = np.min(out1), np.max(out1)
 
     area = coll.dobj[cc][kop]['dgeom']['area']
-    factor = np.ceil((max0 - min0) * (max1 - min1) / area) + 2
+    factor = int(np.ceil((max0 - min0) * (max1 - min1) / area) + 2)
 
     seed0 = np.random.random((nrays * factor,))
     seed1 = np.random.random((nrays * factor,))
@@ -354,7 +411,12 @@ def _seed_optics(
     # imax
     imax = min(nrays*2, iok_pix.size)
 
-    return pts0[iok_pix[:imax]], pts1[iok_pix[:imax]]
+    # optic cent and vectors
+    cent = coll.dobj[cc][kop]['dgeom']['cent']
+    e0 = coll.dobj[cc][kop]['dgeom']['e0']
+    e1 = coll.dobj[cc][kop]['dgeom']['e1']
+
+    return pts0[iok_pix[:imax]], pts1[iok_pix[:imax]], cent, e0, e1
 
 
 # ###############################################################