Repository: ToFuProject/tofu

Issue:
`get_rays_angles_for_single_point_camera2d()` to be implemented

Hints from issue comments:

Solution patch:
diff --git a/tofu/data/_class02_Rays.py b/tofu/data/_class02_Rays.py
index 5ea4b2a9..7147036d 100644
--- a/tofu/data/_class02_Rays.py
+++ b/tofu/data/_class02_Rays.py
@@ -78,6 +78,8 @@ class Rays(Previous):
         key_nseg=None,
         diag=None,
         key_cam=None,
+        # others
+        **kwdargs,
     ):
         """ Add a set of rays
 
@@ -120,6 +122,8 @@ class Rays(Previous):
             key_nseg=key_nseg,
             diag=diag,
             key_cam=key_cam,
+            # kwdargs
+            kwdargs=kwdargs,
         )
 
         # update dicts
@@ -150,20 +154,40 @@ class Rays(Previous):
         self,
         key=None,
         key_cam=None,
+        segment=None,
     ):
-        return _check._get_pts(coll=self, key=key, key_cam=key_cam)
+        """ Return (ptsx, ptsy, ptsz) coordinates for rays
+
+        segment = None => all segments returned, (nseg, ...) array
+        segment = int => (2, ...) array returned
+        """
+        return _check._get_pts(
+            coll=self,
+            key=key,
+            key_cam=key_cam,
+            segment=segment,
+        )
 
     def get_rays_vect(
         self,
         key=None,
         key_cam=None,
         norm=None,
+        segment=None,
     ):
+        """ Return (vx, vy, vz) coordinates for desired rays
+
+        Optionally (default) normalized
+
+        segment = None => all segments returned, (nseg-1, ...) array
+        segment = int => (...,) array returned
+        """
         return _check._get_vect(
             coll=self,
             key=key,
             key_cam=key_cam,
             norm=norm,
+            segment=segment,
         )
 
     # --------------
@@ -376,10 +400,15 @@ class Rays(Previous):
     def add_single_point_camera2d(
         self,
         key=None,
+        # from rays
+        key_rays=None,
+        segment=None,
+        # user-defined
         cent=None,
         nin=None,
         e0=None,
         e1=None,
+        # angles
         angle0=None,
         angle1=None,
         config=None,
@@ -394,19 +423,33 @@ class Rays(Previous):
         """ Add a set of 2d rays from a single point
 
         Rays are sampling a portion of sphere around cent
+
+        cent can be either:
+            - user-defined
+            - derived from intersection of existing rays
+
+        Unit vectors (nin) is also either:
+            - derived from existing rays
+            - user-defined
+
         Portion is defined by 2 angles:
-            - alpha
-            - beta
+            - angle0
+            - angle1
 
         """
 
         return _single_point_cam.main(
             coll=self,
             key=key,
+            # from rays
+            key_rays=key_rays,
+            segment=segment,
+            # user-defined
             cent=cent,
             nin=nin,
             e0=e0,
             e1=e1,
+            # angles
             angle0=angle0,
             angle1=angle1,
             config=config,
@@ -419,6 +462,46 @@ class Rays(Previous):
             units_angles=units_angles,
         )
 
+    # ----------------------
+    # get angles of rasy from single point camera
+    # ----------------------
+
+    def get_rays_angles_from_single_point_camera2d(
+        self,
+        key_single_pt_cam=None,
+        # rays to get angles of
+        key_rays=None,
+        segment=None,
+        # max tolerance
+        tol_radius=None,
+        # optional indices and convex
+        return_indices=None,
+        convex_axis=None,
+        # verb
+        verb=None,
+    ):
+        """ Compute (angle0, angle1) for chosen rays
+
+        As seen from a single point camera2d
+
+        Optionally disable rays with an impact factor too large vs camera
+
+        """
+        return _single_point_cam._get_rays_angles(
+            coll=self,
+            key_single_pt_cam=key_single_pt_cam,
+            # rays to get angles of
+            key_rays=key_rays,
+            segment=segment,
+            # max tolerance
+            tol_radius=tol_radius,
+            # optional indices and convex
+            return_indices=return_indices,
+            convex_axis=convex_axis,
+            # verb
+            verb=verb,
+        )
+
     # ------------------
     # plotting - static
     # ------------------
diff --git a/tofu/data/_class02_single_point_camera.py b/tofu/data/_class02_single_point_camera.py
index ca2c8d66..8872142d 100644
--- a/tofu/data/_class02_single_point_camera.py
+++ b/tofu/data/_class02_single_point_camera.py
@@ -1,6 +1,8 @@
 
 
 import numpy as np
+import scipy.spatial as scpsp
+from matplotlib.path import Path
 import datastock as ds
 
 
@@ -33,10 +35,15 @@ _DANGLES = {
 def main(
     coll=None,
     key=None,
+    # from rays
+    key_rays=None,
+    segment=None,
+    # from user-defined
     cent=None,
     nin=None,
     e0=None,
     e1=None,
+    # angles
     angle0=None,
     angle1=None,
     config=None,
@@ -56,10 +63,15 @@ def main(
     key, cent, nin, e0, e1, dangles, strict = _check(
         coll=coll,
         key=key,
+        # from rays
+        key_rays=key_rays,
+        segment=segment,
+        # user-defined
         cent=cent,
         nin=nin,
         e0=e0,
         e1=e1,
+        # angles
         angle0=angle0,
         angle1=angle1,
         strict=strict,
@@ -86,7 +98,7 @@ def main(
     # angles data
 
     for k0, v0 in dangles.items():
-        if v0['key'] not in coll.dref.keys():
+        if v0['key'] not in coll.ddata.keys():
             coll.add_data(**v0)
 
     # -------------
@@ -134,6 +146,11 @@ def main(
         vect_z=vz,
         # ref
         ref=(dangles['angle0']['ref'], dangles['angle1']['ref']),
+        # specific to ptcam
+        angles=(dangles['angle0']['key'], dangles['angle1']['key']),
+        nin=nin,
+        e0=e0,
+        e1=e1,
         # config
         config=config,
         strict=strict,
@@ -151,6 +168,10 @@ def main(
 def _check(
     coll=None,
     key=None,
+    # from rays
+    key_rays=None,
+    segment=None,
+    # from user-defined
     cent=None,
     nin=None,
     e0=None,
@@ -179,15 +200,71 @@ def _check(
         excluded=lout,
     )
 
+    # ---------------------
+    # rays vs user-defined
+    # ---------------------
+
+    lc = [
+        key_rays is not None,
+        cent is not None,
+    ]
+    if np.sum(lc) != 1:
+        msg = (
+            "Please provide either key_rays xor cent!\n"
+            f"\t- key_rays: {key_rays}\n"
+            f"\t- cent: {cent}\n"
+        )
+        raise Exception(msg)
+
+    # --------------
+    # key_rays
+    # --------------
+
+    if key_rays is not None:
+
+        if isinstance(key_rays, str):
+            key_rays = [key_rays]
+
+        # key_rays
+        lok = list(coll.dobj.get(wrays, {}).keys())
+        key_rays = ds._generic_check._check_var_iter(
+            key_rays, 'key_rays',
+            types=list,
+            types_iter=str,
+            allowed=lok,
+        )
+
+        # segment
+        nseg = [coll.dobj[wrays][kk]['shape'][0] for kk in key_rays]
+        segment = ds._generic_check._check_var(
+            segment, 'segment',
+            types=int,
+            default=0,
+            allowed=list(range(np.min(nseg))),
+        )
+
+        # derive cent
+        cent, nini = _rays_intersection(
+            coll,
+            key_rays=key_rays,
+            segment=segment,
+        )
+
+        # nin
+        if nin is None:
+            nin = nini
+
     # --------------
     # cent
     # --------------
 
-    cent = ds._generic_check._check_flat1darray(
-        cent, 'cent',
-        dtype=float,
-        size=3,
-    )
+    else:
+
+        cent = ds._generic_check._check_flat1darray(
+            cent, 'cent',
+            dtype=float,
+            size=3,
+        )
 
     # --------------
     # unit vectors
@@ -205,7 +282,6 @@ def _check(
     # angle0, angle1 angles
     # --------------
 
-    # angle0
     dangles = {
         'angle0': _check_angles(
             coll=coll,
@@ -333,3 +409,580 @@ def _check_angles(
         'ref': ang_ref,
         'units': ang_units,
     }
+
+
+# ##################################################
+# ##################################################
+#          Get rays intersection
+# ##################################################
+
+
+def _rays_intersection(
+    coll,
+    key_rays=None,
+    segment=None,
+):
+
+    # ----------------
+    # prepare
+    # ----------------
+
+    wrays = coll._which_rays
+    nrays = len(key_rays)
+    nrays_tot = np.array([
+        np.prod(coll.dobj[wrays][kray]['shape'][1:]) for kray in key_rays
+    ])
+
+    cx = np.full((nrays,), np.nan)
+    cy = np.full((nrays,), np.nan)
+    cz = np.full((nrays,), np.nan)
+
+    ux = np.full((nrays,), np.nan)
+    uy = np.full((nrays,), np.nan)
+    uz = np.full((nrays,), np.nan)
+
+    # ----------------
+    # loop on key_rays
+    # ----------------
+
+    for ir, kray in enumerate(key_rays):
+
+        # starting points
+        ptsx, ptsy, ptsz = coll.get_rays_pts(
+            kray,
+            segment=segment,
+        )
+        px = ptsx[0, ...]
+        py = ptsy[0, ...]
+        pz = ptsz[0, ...]
+
+        # unit vectors
+        vx, vy, vz = coll.get_rays_vect(
+            kray,
+            segment=segment,
+        )
+
+        # intersection
+        shape = px.shape + px.shape
+        interx = np.full(shape, np.nan)
+        intery = np.full(shape, np.nan)
+        interz = np.full(shape, np.nan)
+        for i0, ind0 in enumerate(np.ndindex(px.shape)):
+            for i1, ind1 in enumerate(np.ndindex(px.shape)):
+
+                if i1 <= i0:
+                    continue
+
+                sli = ind0 + ind1
+                interx[sli], intery[sli], interz[sli] = _intersect_line2line(
+                    px[ind0], py[ind0], pz[ind0],
+                    px[ind1], py[ind1], pz[ind1],
+                    vx[ind0], vy[ind0], vz[ind0],
+                    vx[ind1], vy[ind1], vz[ind1],
+                )
+
+        cx[ir] = np.nanmean(interx)
+        cy[ir] = np.nanmean(intery)
+        cz[ir] = np.nanmean(interz)
+
+        ux[ir] = np.nanmean(vx)
+        uy[ir] = np.nanmean(vy)
+        uz[ir] = np.nanmean(vz)
+
+    # --------------
+    # take average
+    # --------------
+
+    ntot = np.sum(nrays_tot)
+    cent = np.r_[
+        np.sum(cx * nrays_tot) / ntot,
+        np.sum(cy * nrays_tot) / ntot,
+        np.sum(cz * nrays_tot) / ntot,
+    ]
+    nin = np.r_[
+        np.sum(ux * nrays_tot) / ntot,
+        np.sum(uy * nrays_tot) / ntot,
+        np.sum(uz * nrays_tot) / ntot,
+    ]
+
+    return cent, nin
+
+
+def _intersect_line2line(
+    px0, py0, pz0,
+    px1, py1, pz1,
+    vx0, vy0, vz0,
+    vx1, vy1, vz1,
+):
+
+    # A0M = k0u0
+    # A1N = k1u1
+    # MN = -A0M + A0A1 + A1N
+    # MN = -k0u0 + A0A1 + k1u1
+    # MN.u0 = 0 = -k0 + A0A1.u0 + k1(u0.u1)
+    # MN.u1 = 0 = -k0(u0.u1) + A0A1.u1 + k1
+    #
+    # k0 = A0A1.u0 + k1(u0.u1)
+    # 0 = k1(1 - (u0.u1)^2) - (A0A1.u0)(u0.u1) + A0A1.u1
+    # k1 = ((A0A1.u0)(u0.u1) - A0A1.u1) / (1 - (u0.u1)^2)
+
+    A0A1_x = px1 - px0
+    A0A1_y = py1 - py0
+    A0A1_z = pz1 - pz0
+
+    A0A1_u0 = A0A1_x * vx0 + A0A1_y * vy0 + A0A1_z * vz0
+    A0A1_u1 = A0A1_x * vx1 + A0A1_y * vy1 + A0A1_z * vz1
+    u0_u1 = vx0 * vx1 + vy0 * vy1 + vz0 * vz1
+
+    k1 = ((A0A1_u0) * (u0_u1) - A0A1_u1) / (1. - (u0_u1)**2)
+    k0 = A0A1_u0 + k1 * u0_u1
+
+    # halfway
+    interx = 0.5 * (px0 + k0 * vx0 + px1 + k1 * vx1)
+    intery = 0.5 * (py0 + k0 * vy0 + py1 + k1 * vy1)
+    interz = 0.5 * (pz0 + k0 * vz0 + pz1 + k1 * vz1)
+
+    return interx, intery, interz
+
+
+# ########################################################
+# ########################################################
+#              Get rays angles
+# ########################################################
+
+
+def _get_rays_angles(
+    coll=None,
+    key_single_pt_cam=None,
+    # rays to get angles of
+    key_rays=None,
+    segment=None,
+    # max tolerance
+    tol_radius=None,
+    # optional indices and convex
+    return_indices=None,
+    convex_axis=None,
+    # verb
+    verb=None,
+):
+
+    # -------------
+    # check inputs
+    # -------------
+
+    (
+        key_single_pt_cam, key_rays,
+        segment, tol_radius,
+        return_indices, convex_axis,
+        verb,
+    ) = _check_rays_angles(
+        coll=coll,
+        key_single_pt_cam=key_single_pt_cam,
+        key_rays=key_rays,
+        segment=segment,
+        tol_radius=tol_radius,
+        return_indices=return_indices,
+        convex_axis=convex_axis,
+        verb=verb,
+    )
+
+    # -------------
+    # prepare
+    # -------------
+
+    wrays = coll._which_rays
+    ref = coll.dobj[wrays][key_rays]['ref'][1:]
+    shape = coll.dobj[wrays][key_rays]['shape'][1:]
+    kang0, kang1 = coll.dobj[wrays][key_single_pt_cam]['angles']
+    units = coll.ddata[kang0]['units']
+
+    # unit vectors
+    nin = coll.dobj[wrays][key_single_pt_cam]['nin']
+    e0 = coll.dobj[wrays][key_single_pt_cam]['e0']
+    e1 = coll.dobj[wrays][key_single_pt_cam]['e1']
+
+    ptsx, ptsy, ptsz = coll.get_rays_pts(key_single_pt_cam)
+    cent = np.r_[ptsx[0, 0, 0], ptsy[0, 0, 0], ptsz[0, 0, 0]]
+
+    # -------------
+    # get impact factor vs tol_radius
+    # -------------
+
+    ptsx, ptsy, ptsz = coll.get_rays_pts(key_rays, segment=segment)
+    vx, vy, vz = coll.get_rays_vect(key_rays, segment=segment, norm=True)
+
+    p0c_x = (ptsx[0, ...] - cent[0])
+    p0c_y = (ptsy[0, ...] - cent[1])
+    p0c_z = (ptsz[0, ...] - cent[2])
+
+    crossx = p0c_y * vz - p0c_z * vy
+    crossy = p0c_z * vx - p0c_x * vz
+    crossz = p0c_x * vy - p0c_y * vx
+
+    impact = np.sqrt(crossx**2 + crossy**2 + crossz**2)
+
+    # verb
+    if verb is True:
+        msg = (
+            "Angles computation for rays, seen from single_point_camera2d:\n"
+            f"\t- single_point_camera2d: {key_single_pt_cam}\n"
+            f"\t- key_rays: {key_rays}\n"
+            f"\t- max impact: {np.max(impact)} m\n"
+        )
+        print(msg)
+
+    # -------------
+    # compute
+    # -------------
+
+    # intialize
+    ang0 = np.full(shape, np.nan)
+    ang1 = np.full(shape, np.nan)
+
+    # vect from cent
+    vx = ptsx[1, ...] - cent[0]
+    vy = ptsy[1, ...] - cent[1]
+    vz = ptsz[1, ...] - cent[2]
+    vnorm = np.sqrt(vx**2 + vy**2 + vz**2)
+    vx = vx / vnorm
+    vy = vy / vnorm
+    vz = vz / vnorm
+
+    # derive angles
+    sin1 = vx * e1[0] + vy * e1[1] + vz * e1[2]
+    cos1_sin0 = vx * e0[0] + vy * e0[1] + vz * e0[2]
+    cos1_cos0 = vx * nin[0] + vy * nin[1] + vz * nin[2]
+    cos12 = cos1_sin0**2 + cos1_cos0**2
+    cos1 = np.sqrt(cos12)
+    ang1 = np.arctan2(sin1, cos1)
+    ang0 = np.arctan2(cos1_sin0 / cos1, cos1_cos0 / cos1)
+
+    # adjust units
+    if str(units) == 'deg':
+        ang0 = ang0 * 180/np.pi
+        ang1 = ang1 * 180/np.pi
+
+    # -------------
+    # indices
+    # -------------
+
+    if return_indices is True:
+
+        bin0 = coll.ddata[kang0]['data']
+        dang0 = bin0[1] - bin0[0]
+        bin0 = 0.5*(bin0[1:] + bin0[:-1])
+        bin0 = np.r_[bin0[0]-dang0, bin0]
+
+        bin1 = coll.ddata[kang0]['data']
+        dang1 = bin1[1] - bin1[0]
+        bin1 = 0.5*(bin1[1:] + bin1[:-1])
+        bin1 = np.r_[bin1[0]-dang1, bin1]
+
+        ind0 = np.searchsorted(bin0, ang0) - 1
+        ind1 = np.searchsorted(bin1, ang1) - 1
+
+        # adjust for edges
+        ind0[ind0 < 0] = 0
+        ind1[ind1 < 0] = 0
+
+        # -------
+        # convex_axis
+
+        if convex_axis is not False:
+
+            # shape_angles
+            shape_angles = coll.dobj[wrays][key_single_pt_cam]['shape'][1:]
+
+            # shape
+            shape_c = tuple([
+                ss for ii, ss in enumerate(shape)
+                if ii not in convex_axis
+            ])
+
+            # unique call
+            if len(shape_c) == 0:
+                sli = tuple([slice(None) for ii in shape])
+                dhull = {
+                    (0,): _convexhull(
+                        shape_angles=shape_angles,
+                        ind0=ind0,
+                        ind1=ind1,
+                        sli=sli,
+                    )
+                }
+
+            else:
+                # initialize
+                dhull = {}
+                sli = np.array([
+                    slice(None) if ii in convex_axis else 0
+                    for ii in range(len(shape))
+                ])
+                isli = np.array([
+                    ii for ii in range(len(shape))
+                    if ii not in convex_axis
+                ])
+
+                i0 = np.arange(shape_angles[0])
+                i1 = np.arange(shape_angles[1])
+                i01_full = np.array([
+                    np.repeat(i0[:, None], i1.size, axis=1).ravel(),
+                    np.repeat(i1[None, :], i0.size, axis=0).ravel(),
+                ]).T
+
+                # loop
+                for ii, ind in enumerate(np.ndindex(shape_c)):
+
+                    # update slice
+                    sli[isli] = ind
+                    slii = tuple(sli)
+
+                    # get convex_axis hull
+                    dhull[ind] = _convexhull(
+                        shape_angles=shape_angles,
+                        ind0=ind0,
+                        ind1=ind1,
+                        sli=slii,
+                        i01_full=i01_full,
+                    )
+
+    # -------------
+    # output
+    # -------------
+
+    dout = {
+        'angle0': {
+            'key': f'{key_rays}_{key_single_pt_cam}_ang0',
+            'data': ang0,
+            'units': units,
+            'ref': ref,
+            'dim': 'angle',
+        },
+        'angle1': {
+            'key': f'{key_rays}_{key_single_pt_cam}_ang1',
+            'data': ang1,
+            'units': units,
+            'ref': ref,
+            'dim': 'angle',
+        },
+        'impact': {
+            'key': f'{key_rays}_{key_single_pt_cam}_impact',
+            'data': impact,
+            'ref': ref,
+            'units': 'm',
+            'dim': 'distance',
+        },
+        'indices': {},
+    }
+
+    if return_indices:
+        dout['ind0'] = {
+            'key': f'{key_rays}_{key_single_pt_cam}_ind0',
+            'data': ind0,
+            'ref': ref,
+        }
+        dout['ind1'] = {
+            'key': f'{key_rays}_{key_single_pt_cam}_ind1',
+            'data': ind1,
+            'ref': ref,
+        }
+
+        if convex_axis is not False:
+            dout['hull'] = dhull
+
+    return dout
+
+
+# ########################################################
+# ########################################################
+#              check inputs rays angles
+# ########################################################
+
+
+def _check_rays_angles(
+    coll=None,
+    key_single_pt_cam=None,
+    # rays to get angles of
+    key_rays=None,
+    segment=None,
+    # max tolerance
+    tol_radius=None,
+    # optional indices
+    return_indices=None,
+    convex_axis=None,
+    # verb
+    verb=None,
+):
+
+    # -------------
+    # keys
+    # -------------
+
+    wrays = coll._which_rays
+
+    # key_rays
+    lok = list(coll.dobj.get(wrays, {}).keys())
+    key_rays = ds._generic_check._check_var(
+        key_rays, 'key_rays',
+        types=str,
+        allowed=lok,
+    )
+
+    # key_single_pt_cam
+    lok = [
+        k0 for k0 in lok
+        if isinstance(coll.dobj[wrays][k0].get('angles'), tuple)
+        and len(coll.dobj[wrays][k0]['angles']) == 2
+        and all([
+            ss in coll.ddata.keys()
+            for ss in coll.dobj[wrays][k0].get('angles')
+        ])
+        and isinstance(coll.dobj[wrays][k0].get('nin'), np.ndarray)
+        and isinstance(coll.dobj[wrays][k0].get('e0'), np.ndarray)
+        and isinstance(coll.dobj[wrays][k0].get('e1'), np.ndarray)
+    ]
+    key_single_pt_cam = ds._generic_check._check_var(
+        key_single_pt_cam, 'key_single_pt_cam',
+        types=str,
+        allowed=lok,
+    )
+
+    # -----------------
+    # segment
+    # -----------------
+
+    segment = int(ds._generic_check._check_var(
+        segment, 'segment',
+        types=(int, float),
+        default=0,
+    ))
+
+    # -----------------
+    # tol_radius
+    # -----------------
+
+    if tol_radius is not None:
+        tol_radius = float(ds._generic_check._check_var(
+            tol_radius, 'tol_radius',
+            types=(int, float),
+            sign='>0',
+        ))
+
+    # -----------------
+    # return_indices
+    # -----------------
+
+    return_indices = ds._generic_check._check_var(
+        return_indices, 'return_indices',
+        types=bool,
+        default=convex_axis not in [None, False],
+    )
+
+    # --------------
+    # convex_axis
+    # --------------
+
+    if np.isscalar(convex_axis) and not isinstance(convex_axis, bool):
+        convex_axis = (int(convex_axis),)
+
+    convex_axis = ds._generic_check._check_var(
+        convex_axis, 'convex_axis',
+        types=(tuple, bool),
+        default=False,
+    )
+
+    # case by case
+    ndim = len(coll.dobj[wrays][key_rays]['shape'][1:])
+    if convex_axis is True:
+        convex_axis = tuple(range(ndim))
+
+    elif convex_axis is not False:
+        convex_axis = ds._generic_check._check_var_iter(
+            convex_axis, 'convex_axis',
+            types=tuple,
+            types_iter=int,
+            allowed=list(range(ndim)) + list(range(-1, -ndim-1, -1)),
+        )
+
+        convex_axis = np.array(convex_axis)
+        ineg = convex_axis < 0
+        convex_axis[ineg] = ndim + convex_axis[ineg]
+        convex_axis = tuple(sorted(convex_axis))
+
+    # -----------------
+    # verb
+    # -----------------
+
+    verb = ds._generic_check._check_var(
+        verb, 'verb',
+        types=bool,
+        default=True,
+    )
+
+    return (
+        key_single_pt_cam, key_rays,
+        segment, tol_radius,
+        return_indices, convex_axis,
+        verb,
+    )
+
+
+# ####################################################
+# ####################################################
+#            get ConvexHull
+# ####################################################
+
+
+def _convexhull(
+    shape_angles=None,
+    ind0=None,
+    ind1=None,
+    sli=None,
+    i01_full=None,
+):
+    # ---------------
+    # initialize
+    # ---------------
+
+    ind = np.zeros(shape_angles, dtype=bool)
+
+    # ---------------
+    # get convex hull
+    # ---------------
+
+    ui0 = np.unique(ind0[sli])
+    ui1 = np.unique(ind1[sli])
+    if ui0.size == 1:
+        npts = ui1[-1] - ui1[0] + 1
+        slii = (
+            np.full((npts,), ui0[0]),
+            np.arange(ui1[0], min(ui1[-1]+1, ind.shape[1])),
+        )
+        ind[slii] = True
+
+    elif ui1.size == 1:
+        npts = ui0[-1] - ui0[0] + 1
+        slii = (
+            np.arange(ui0[0], min(ui0[-1]+1, ind.shape[0])),
+            np.full((npts,), ui1[0]),
+        )
+        ind[slii] = True
+
+    # ---------------
+    # get convex hull
+    # ---------------
+
+    else:
+        i0_flat = ind0[sli].ravel()
+        i1_flat = ind1[sli].ravel()
+        hull = scpsp.ConvexHull(
+            np.array([i0_flat, i1_flat]).T
+        )
+
+        # get indices of angles in hull
+        path = Path(
+            np.array([i0_flat[hull.vertices], i1_flat[hull.vertices]]).T
+        )
+        ind = path.contains_points(i01_full).reshape(shape_angles)
+
+    return ind
diff --git a/tofu/data/_class08_generate_rays.py b/tofu/data/_class08_generate_rays.py
index 8b2aba43..038f03e0 100644
--- a/tofu/data/_class08_generate_rays.py
+++ b/tofu/data/_class08_generate_rays.py
@@ -204,6 +204,7 @@ def _check(
                 f"Arg '{kdict}' must be a dict with at least one of:\n"
                 + "\n".join(lstr)
                 + f"\nFed to {kfunc}\n"
+                + f"\nProvided:\n{vdict}\n"
             )
             raise Exception(msg)
 
@@ -267,7 +268,7 @@ def _check(
 
         else:
             if len(key_cam) == 1 and isinstance(key_rays, str):
-                key_rays[key_cam[0]] = key_rays
+                key_rays = {key_cam[0]: key_rays}
 
         ncam = len(key_cam)
         c0 = (
diff --git a/tofu/data/_class2_check.py b/tofu/data/_class2_check.py
index aaca015c..edb5e6db 100644
--- a/tofu/data/_class2_check.py
+++ b/tofu/data/_class2_check.py
@@ -8,10 +8,10 @@ import datastock as ds
 from ..geom import CamLOS1D
 
 
-# ################################################################
-# ################################################################
+# ############################################################
+# ############################################################
 #                   Rays
-# ################################################################
+# ############################################################
 
 
 def _start_vect(ax=None, ay=None, az=None, name=None):
@@ -278,8 +278,9 @@ def _check_inputs(
         )
 
         # lspectro
+        loptics = coll.dobj['diagnostic'][diag]['doptics'][key_cam]['optics']
         lspectro = [
-            oo for oo in coll.dobj['diagnostic'][diag]['doptics'][key_cam]['optics']
+            oo for oo in loptics
             if oo in coll.dobj.get('crystal', {}).keys()
             or oo in coll.dobj.get('grating', {}).keys()
         ]
@@ -297,10 +298,10 @@ def _check_inputs(
     )
 
 
-# ################################################################
-# ################################################################
+# ############################################################
+# ############################################################
 #                   Main
-# ################################################################
+# ############################################################
 
 
 def _rays(
@@ -335,6 +336,8 @@ def _rays(
     key_nseg=None,
     diag=None,
     key_cam=None,
+    # kwdargs
+    kwdargs=None,
 ):
 
     # -------------
@@ -404,12 +407,6 @@ def _rays(
             vy = vy / norm
             vz = vz / norm
 
-            # sca = (
-                # vx[1:, ...] * vx[:-1, ...]
-                # + vy[1:, ...] * vy[:-1, ...]
-                # + vz[1:, ...] * vz[:-1, ...]
-            # )
-
     # -----------------------------
     # compute ray-tracing
     # -----------------------------
@@ -600,13 +597,14 @@ def _rays(
         key_nseg=key_nseg,
         ref=ref,
         shaperef=shaperef,
+        kwdargs=kwdargs,
     )
 
 
-# ##################################################################
-# ##################################################################
+# ###########################################################
+# ###########################################################
 #                   Rays - check key
-# ##################################################################
+# ###########################################################
 
 
 def _check_key(coll=None, key=None, key_cam=None):
@@ -646,10 +644,38 @@ def _check_key(coll=None, key=None, key_cam=None):
     return kray
 
 
-# ##################################################################
-# ##################################################################
+# #########################################################
+# #########################################################
+#                   Rays - check segment
+# #########################################################
+
+
+def _check_segment(coll, key, segment):
+
+    if segment is not None:
+
+        wrays = coll._which_rays
+        nseg = coll.dobj[wrays][key]['shape'][0]
+        lok = list(range(nseg)) + [-1]
+        segment = int(ds._generic_check._check_var(
+            segment, 'segment',
+            types=(int, float),
+        ))
+        if segment not in lok:
+            msg = (
+                f"Arg segment for rays '{key}' must be in:\n"
+                f"\t- allowed: {lok}\n"
+                f"\t- Provided: {segment}\n"
+            )
+            raise Exception(msg)
+
+    return segment
+
+
+# #########################################################
+# #########################################################
 #                   Rays - get start and vect
-# ##################################################################
+# #########################################################
 
 
 def _get_start(
@@ -680,12 +706,14 @@ def _get_pts(
     coll=None,
     key=None,
     key_cam=None,
+    segment=None,
 ):
 
     # ---------
     # check key
 
     key = _check_key(coll=coll, key=key, key_cam=key_cam)
+    segment = _check_segment(coll, key, segment)
 
     # ---------
     # get start
@@ -702,9 +730,20 @@ def _get_pts(
     ptsz = coll.ddata[ptsz]['data']
 
     # concatenate
-    ptsx = np.insert(ptsx, 0, stx, axis=0)
-    ptsy = np.insert(ptsy, 0, sty, axis=0)
-    ptsz = np.insert(ptsz, 0, stz, axis=0)
+    if segment is None:
+        ptsx = np.insert(ptsx, 0, stx, axis=0)
+        ptsy = np.insert(ptsy, 0, sty, axis=0)
+        ptsz = np.insert(ptsz, 0, stz, axis=0)
+
+    elif segment == 0:
+        ptsx = np.insert(ptsx[0:1, ...], 0, stx, axis=0)
+        ptsy = np.insert(ptsy[0:1, ...], 0, sty, axis=0)
+        ptsz = np.insert(ptsz[0:1, ...], 0, stz, axis=0)
+
+    else:
+        ptsx = ptsx[segment-1:segment+1, ...]
+        ptsy = ptsy[segment-1:segment+1, ...]
+        ptsz = ptsz[segment-1:segment+1, ...]
 
     return ptsx, ptsy, ptsz
 
@@ -714,12 +753,14 @@ def _get_vect(
     key=None,
     key_cam=None,
     norm=None,
+    segment=None,
 ):
 
     # ---------
     # check key
 
     key = _check_key(coll=coll, key=key, key_cam=key_cam)
+    segment = _check_segment(coll, key, segment)
 
     # norm
     norm = ds._generic_check._check_var(
@@ -754,13 +795,19 @@ def _get_vect(
         vy = vy / norm
         vz = vz / norm
 
+    # segment
+    if segment is not None:
+        vx = vx[segment, ...]
+        vy = vy[segment, ...]
+        vz = vz[segment, ...]
+
     return vx, vy, vz
 
 
-# ##################################################################
-# ##################################################################
+# ##########################################################
+# ##########################################################
 #                   store
-# ##################################################################
+# ##########################################################
 
 
 def _make_dict(
@@ -782,6 +829,8 @@ def _make_dict(
     key_nseg=None,
     ref=None,
     shaperef=None,
+    # kwdargs
+    kwdargs=None,
 ):
 
     # --------
@@ -801,7 +850,6 @@ def _make_dict(
         )
         raise Exception(msg)
 
-
     assert nseg == reflections_nb + 1 + nextra
 
     # key_nseg
@@ -993,6 +1041,7 @@ def _make_dict(
                 'alpha': kalpha,
                 'reflect_dalpha': kdalpha,
                 'reflect_dbeta': kdbeta,
+                **kwdargs,
             },
         },
     }


Test patch:
diff --git a/tofu/tests/tests08_diagnostics/test_01_diagnostics.py b/tofu/tests/tests08_diagnostics/test_01_diagnostics.py
index 1feaa02d..8d37dacf 100644
--- a/tofu/tests/tests08_diagnostics/test_01_diagnostics.py
+++ b/tofu/tests/tests08_diagnostics/test_01_diagnostics.py
@@ -204,7 +204,7 @@ def _apertures():
     # kl0 = 0.02 * np.linspace(-1, 1, npix)
     klin = 0.2
     delta = 0.05
-    cents_x = cent[0] - (klin - delta) * nins[:, 0]    #  + kl0 * e0[0]
+    cents_x = cent[0] - (klin - delta) * nins[:, 0]    # + kl0 * e0[0]
     cents_y = cent[1] - (klin - delta) * nins[:, 1]    # + kl0 * e0[1]
     cents_z = cent[2] - (klin - delta) * nins[:, 2]    # + kl0 * e0[2]
 
@@ -542,7 +542,9 @@ def _diagnostics():
 
     # d12: 1d collimator camera, one aperture per pixel
     d12 = {
-        'doptics': tuple(['cam0000'] + [k1 for k1 in [f'lap{ii}' for ii in range(10)]]),
+        'doptics': tuple(
+            ['cam0000'] + [k1 for k1 in [f'lap{ii}' for ii in range(10)]]
+        ),
     }
 
     # d13: 1d collimator-hybrid camera
@@ -550,7 +552,10 @@ def _diagnostics():
         'doptics': {
             'cam00000': {
                 'optics': [f'lap{ii}' for ii in range(10)] + ['ap0'],
-                'paths': np.concatenate((np.eye(10), np.ones((10, 1))), axis=1),
+                'paths': np.concatenate(
+                    (np.eye(10), np.ones((10, 1))),
+                    axis=1,
+                ),
             },
         },
     }
@@ -742,7 +747,10 @@ class Test01_Diagnostic():
         lrays = list(self.coll.dobj['rays'].keys())
         ldiag = [
             k0 for k0, v0 in self.coll.dobj['diagnostic'].items()
-            if any([v1.get('los') is not None for v1 in v0['doptics'].values()])
+            if any([
+                v1.get('los') is not None
+                for v1 in v0['doptics'].values()
+            ])
         ]
         lk = [lrays] + ldiag
         for ii, k0 in enumerate(lk):
@@ -816,7 +824,7 @@ class Test01_Diagnostic():
             if self.coll.dobj['diagnostic'][k0]['spectro']:
                 continue
 
-            if self.coll.dobj['diagnostic'][k0]['doptics'][lcam[0]]['dvos'].get('keym') is None:
+            if doptics[lcam[0]]['dvos'].get('keym') is None:
                 self.coll.compute_diagnostic_vos(
                     # keys
                     key_diag=k0,
@@ -836,14 +844,53 @@ class Test01_Diagnostic():
         plt.close('all')
 
     def test07_add_single_point_camera2d(self):
+
         self.coll.add_single_point_camera2d(
-            cent=np.r_[3, 0, 0],
-            nin=np.r_[-1, 0, 1],
-            angle0=25,
-            angle1=25,
+            key='ptcam',
+            key_rays='diag5_cam22_los',
+            angle0=55,
+            angle1=55,
             config=self.conf_touch,
         )
 
+        # add rays
+        dsamp = {'dedge': {'res': 'min'}, 'dsurface': {'nb': 3}}
+        self.coll.add_rays_from_diagnostic(
+            key='diag5',
+            dsampling_pixel=dsamp,
+            dsampling_optics=dsamp,
+            optics=-1,
+            config=self.conf_touch,
+            store=True,
+            strict=None,
+            key_rays='diag5_cam22_rays',
+            overwrite=None,
+        )
+
+        # get angles from rays
+        dout = self.coll.get_rays_angles_from_single_point_camera2d(
+            key_single_pt_cam='ptcam',
+            key_rays='diag5_cam22_los',
+            return_indices=False,
+        )
+
+        # get angles from rays
+        dout = self.coll.get_rays_angles_from_single_point_camera2d(
+            key_single_pt_cam='ptcam',
+            key_rays='diag5_cam22_rays',
+            return_indices=True,
+            convex_axis=False,
+        )
+
+        # get angles from rays
+        dout = self.coll.get_rays_angles_from_single_point_camera2d(
+            key_single_pt_cam='ptcam',
+            key_rays='diag5_cam22_rays',
+            return_indices=True,
+            convex_axis=(-1, -2),
+        )
+        assert isinstance(dout, dict)
+
     def test07_add_rays_from_diagnostic(self):
         for ii, (k0, v0) in enumerate(self.coll.dobj['diagnostic'].items()):
             noptics = any([