Implemented method for finding which Bravais lattice points lie withi…

…n a bounding box in n-dimensions

polymerist/maths/lattices/coordops.py

@@ -6,6 +6,7 @@
 import numpy as np
 from itertools import product as cartesian_product
 
+from .integer import CubicIntegerLattice
 from ...genutils.typetools.numpytypes import Shape, N, M, DType
 from ...genutils.typetools.categorical import Numeric
 T = TypeVar('T', bound=Numeric)
@@ -99,13 +100,41 @@ def nearest_int_coord_along_normal(point : np.ndarray[Shape[N], Numeric], normal
         int_point = point 
     else:
         int_bbox = BoundingBox(np.vstack([min_int_bound_point, max_int_bound_point]))  # generate bounding box from extremal positions
-        dots = np.inner((int_bbox.extrema - point), normal) # take dot product between the normal and the direction vectors from the control point to the integer bounding points
+        dots = np.inner((int_bbox.vertices - 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
-        furthest_point, similarity = int_bbox.extrema[i], dots[i]
+        furthest_point, similarity = int_bbox.vertices[i], dots[i]
 
         if similarity <= 0:
             raise ValueError(f'Could not locate valid integral point in normal direction (closest found was {furthest_point} with dot product {similarity})')
         else:
             int_point = furthest_point 
 
-    return int_point.astype(int)
+    return int_point.astype(int)
+
+def identify_bravais_points_within_bbox(lattice_vectors : np.ndarray[Shape[N, N], T], bbox : BoundingBox) -> tuple[np.ndarray[Shape[M, N], T], CubicIntegerLattice]:
+    '''Locate all lattice points generated by a set of Bravais lattice vectors in N-dimensions which fall within a given bounding box
+    Returns the coordinate vector of the internal lattice points and a CubicIntegerLattice containing the corresponding lattice vector multiplicities'''
+    # 1) transform the bounding box into the inverse domain, where Bravais lattice points "look like" integer points
+    inv_bbox_vertices = bbox.vertices @ np.linalg.inv(lattice_vectors.T)
+    inv_bbox_center = COM(inv_bbox_vertices)
+    directors = (inv_bbox_vertices - inv_bbox_center) # directional vectors at each vertex which are away from transformed box center
+    directors /= np.linalg.norm(directors, axis=1)[:, None] # normalize direction vectors
+
+    # 2) find the nearest integer-valued points to transformed bounding box points in directions "away" from the bounding box
+    integral_bounding_points = np.vstack([ 
+        nearest_int_coord_along_normal(point, director) # ...(ensures no bounded points end up outside of final box)
+            for point, director in zip(inv_bbox_vertices, directors)
+    ])
+
+    # 3) produce collection of integer point which contains AT LEAST the integral points bounded by the transformed nearest-integrer bounding points
+    trial_bbox = BoundingBox(integral_bounding_points)  # TODO: box is pretty inefficient for highly-sheared unit vectors: need simple method to test for integral points within a bounding volume (or simplex) 
+    trial_int_latt = CubicIntegerLattice(trial_bbox.dims + 1) # need K+1 unit-spaced points to space an interval K units long (e.g. [0-2] -> (0, 1, 2))
+    trial_int_latt.points += trial_bbox.minimum # offset lattice from default at origin
+
+    # 4) transforms the trial points back to the normal domain and performa final enclosure check
+    trial_latt_points = trial_int_latt.linear_transformation(lattice_vectors)
+    is_inside = bbox.surrounds(trial_latt_points).all(axis=1) # mask of which points have ALL coordinates inside the bounding box
+    inside_latt_points    = trial_latt_points[is_inside]
+    trial_int_latt.points = trial_int_latt.points[is_inside]  # screen out integral points which fall outside the bounding box when mapped to lattice points
+
+    return inside_latt_points, trial_int_latt