From ce8aaf696cbf360a3ed63c56f1ce90dd9ea7f3a9 Mon Sep 17 00:00:00 2001
From: Timotej Bernat <tibe3324@colorado.edu>
Date: Wed, 24 Apr 2024 20:02:37 -0600
Subject: [PATCH] Implemented calculation for nearest integer-valued point in a
 given normal direction relative to a control point

---
 polymerist/maths/lattices/coordops.py | 36 +++++++++++++++++++++++++--
 1 file changed, 34 insertions(+), 2 deletions(-)

diff --git a/polymerist/maths/lattices/coordops.py b/polymerist/maths/lattices/coordops.py
index f313be3..18833a0 100644
--- a/polymerist/maths/lattices/coordops.py
+++ b/polymerist/maths/lattices/coordops.py
@@ -7,8 +7,10 @@
 from itertools import product as cartesian_product
 
 from ...genutils.typetools.numpytypes import Shape, N, M, DType
+from ...genutils.typetools.categorical import Numeric
 
 
+# COORDINATE-TO-POINT DISTANCES
 def mean_coord(coords : np.ndarray[Shape[M, N], DType], weights : Optional[np.ndarray[Shape[M, 1], DType]]=None) -> np.ndarray[Shape[N], DType]:
     '''The average (center-of-mass) coordinate of a vector of coordinates, with optional array of weightsfor each coordinates'''
     if weights is None:
@@ -27,7 +29,7 @@ def dists_to_centroid(coords : np.ndarray[Shape[M, N], DType], norm_order : Opti
     '''The distance of each coordinate in an array of coordinates to the coordinates' centroid'''
     return dists_to_point(coords, point=mean_coord(coords), norm_order=norm_order)
 
-
+# BOUNDING BOXES
 def bounding_box_dims(coords : np.ndarray[Shape[M, N], DType]) -> np.ndarray[Shape[N], DType]:
     '''The N-dimensional tight bounding box side lengths around an array of coordinates'''
     return coords.ptp(axis=0)
@@ -47,4 +49,34 @@ def coords_inside_bbox(coords : np.ndarray[Shape[M, N], DType], lower : np.ndarr
     '''Boolean mask of whether coordinates are within the boundaries of some bounding box
     With strict=True, points on the boundary are not considered inside; with strict=False, points on the boundary are considered insde as well'''
     less_funct = np.less if strict else np.less_equal # set "<" vs "<=" check by strict flag
-    return less_funct(lower, coords) & less_funct(coords, upper)
\ No newline at end of file
+    return less_funct(lower, coords) & less_funct(coords, upper)
+
+# NORMAL VECTORS
+def nearest_int_coord_along_normal(point : np.ndarray[Shape[N], Numeric], normal : np.ndarray[Shape[N], Numeric]) -> np.ndarray[Shape[N], int]:
+    '''
+    Takes an N-dimensional control point and an N-dimensional normal vector
+    Returns the integer-valued point nearest to the control point which lies in
+    the normal direction relative to the control point
+    '''
+    # Check that inputs have vector-like shapes and compatible dimensions
+    assert(point.ndim == 1)
+    assert(normal.ndim == 1)
+    assert(point.size == normal.size)
+
+    min_int_bound_point = np.floor(point)
+    max_int_bound_point = np.ceil( point)
+    if np.isclose(min_int_bound_point, max_int_bound_point).all(): # edge case: if already on an integer-valued point, the fllor and ceiling will be identical
+        int_point = point 
+    else:
+        int_bound_extrema = np.vstack([min_int_bound_point, max_int_bound_point]) # package extremal pointsinto single array
+        int_bound_points = bounding_box_points(int_bound_extrema) # expand extremal points into complete N-dimensional integer box around the control point
+        dots = np.inner((int_bound_points - point), normal) # take dot product between the normal and the direction vectors from the control point to the integer bounding points
+        i = dots.argmax() # position of integer point in most similar direction to normal
+
+        if dots[i] > 0:
+            int_point = int_bound_points[i]
+        else:
+            raise ValueError(f'Could not locate valid integral point in normal direction (closest found was {int_bound_points[i]} with dot product {dots[i]})')
+        pass
+
+    return int_point.astype(int)
\ No newline at end of file