wip: xform problems box2box

bdm_voxel_builder/grid/base.py

@@ -8,11 +8,17 @@
 import numpy.typing as npt
 import pyopenvdb as vdb
 from compas.colors import Color
+from scipy.spatial import QhullError
 
 from bdm_voxel_builder import TEMP_DIR
-from bdm_voxel_builder.helpers.numpy import convert_array_to_pts
-from bdm_voxel_builder.helpers.savepaths import get_savepath
-from bdm_voxel_builder.helpers.vdb import xform_to_compas, xform_to_vdb
+from bdm_voxel_builder.helpers import (
+    convert_array_to_pts,
+    get_linear_xform_between_2_boxes,
+    get_savepath,
+    xform_to_compas,
+    xform_to_vdb,
+)
+from bdm_voxel_builder.helpers.pointcloud import pointcloud_to_grid_array
 
 
 class Grid:
@@ -61,7 +67,8 @@ def grid_size(self, value):
 
     def get_local_bbox(self) -> cg.Box:
         """Returns a bounding box containing the grid, 0, 0, 0 to ijk"""
-        return cg.Box.from_diagonal(((0, 0, 0), self.grid_size))
+        xsize, ysize, zsize = self.grid_size
+        return cg.Box.from_diagonal(cg.Line((0, 0, 0), (xsize, ysize, zsize)))
 
     def get_world_bbox(self) -> cg.Box:
         return self.get_local_bbox().transformed(self.xform)
@@ -102,7 +109,7 @@ def get_value_at_index(self, index=(0, 0, 0)):
         i, j, k = index
         return self.array[i][j][k]
 
-    def get_active_voxels(self, array):
+    def get_active_voxels(self):
         """returns indicies of nonzero values
         list of coordinates
             shape = [3,n]"""
@@ -163,34 +170,33 @@ def from_vdb(cls, grid: os.PathLike | vdb.GridBase, name: str = None):
         )
 
     @classmethod
-    def from_pointcloud(cls, pointcloud: cg.Pointcloud, name: str = None):
-        vdb.grid.createLevelSetFromPolygons(pointcloud.points)
-        if isinstance(grid, os.PathLike):
-            grids = vdb.readAllGridMetadata(str(grid))
-
-            if not name and len(grids) > 1:
-                print(
-                    "File contains more than one grid, ",
-                    f"only processing first named {grids[0].name}",
-                )
+    def from_pointcloud(
+        cls,
+        pointcloud: cg.Pointcloud,
+        grid_size: list[int, int, int] | int,
+        name: str = None,
+    ):
+        try:
+            bbox_pointcloud = pointcloud.obb
+        except QhullError:
+            bbox_pointcloud = pointcloud.aabb
 
-            name = name or grids[0].name
-            grid = vdb.read(str(grid), name)
+
 
-        bbox_min = grid.metadata["file_bbox_min"]
-        bbox_max = grid.metadata["file_bbox_max"]
+        # -1 because the grid is 0 indexed
+        if isinstance(grid_size, int):
+            isize = jsize = ksize = grid_size - 1
+        else:
+            isize, jsize, ksize = (n - 1 for n in grid_size)
 
-        shape = np.array(bbox_max) - np.array(bbox_min)
-        arr = np.zeros(shape)
+        bbox_grid_frame = cg.Frame(point=(isize / 2, jsize / 2, ksize / 2))
 
-        # rotate the grid to make Z up
-        grid.transform.rotate(math.pi / 2, vdb.Axis.X)
+        bbox_grid = cg.Box(xsize=isize, ysize=jsize, zsize=ksize, frame=bbox_grid_frame)
 
-        grid.copyToArray(arr, ijk=bbox_min)
+        xform = get_linear_xform_between_2_boxes(bbox_pointcloud, bbox_grid)
 
-        return cls(
-            grid_size=arr.shape,
-            name=name or grid.name,
-            array=arr,
-            xform=xform_to_compas(grid.transform),
+        array = pointcloud_to_grid_array(
+            pointcloud=pointcloud.transformed(xform), grid_size=grid_size
         )
+
+        return cls(grid_size=grid_size, name=name, xform=xform, array=array)

bdm_voxel_builder/helpers/__init__.py

@@ -0,0 +1,31 @@
+# ruff: noqa: F401, F403
+
+from .common import get_nth_newest_file_in_folder
+from .compas import box_from_corner_frame, get_linear_xform_between_2_boxes
+from .math import calculate_chance, remap, remap_trim
+from .numpy import (
+    NB_INDEX_DICT,
+    clip_indices_to_grid_size,
+    conditional_fill,
+    convert_array_to_pts,
+    create_random_array,
+    crop_array,
+    get_cube_array_indices,
+    get_mask_zone_xxyyzz,
+    get_sub_array,
+    make_solid_box_xxyyzz,
+    make_solid_box_xxz,
+    make_solid_box_z,
+    random_choice_index_from_best_n,
+    random_choice_index_from_best_n_old,
+    save_ndarray,
+    sort_pts_by_values,
+)
+from .pointcloud import (
+    ply_to_compas,
+    ply_to_numpy,
+    pointcloud_from_ndarray,
+    save_pointcloud,
+)
+from .savepaths import get_savepath
+from .vdb import xform_to_compas, xform_to_vdb

bdm_voxel_builder/helpers/compas.py

@@ -0,0 +1,33 @@
+import compas.geometry as cg
+
+
+def box_from_corner_frame(frame: cg.Frame, xsize:float, ysize:float, zsize: float):
+    """Create a box from origin and size."""
+    center_pt = frame.point.copy()
+    center_pt += cg.Vector(xsize / 2, ysize / 2, zsize / 2)
+
+    center_frame = cg.Frame(center_pt, xaxis=frame.xaxis, yaxis=frame.yaxis)
+
+    return cg.Box(xsize=xsize, ysize=ysize, zsize=zsize, frame=center_frame)    
+
+
+def get_linear_xform_between_2_boxes(from_box: cg.Box, to_box: cg.Box):
+    """Get the linear transformation between two bounding boxes."""
+    f_pt = from_box.frame.point.copy()
+    f_pt -= cg.Point(from_box.width / 2, from_box.height / 2, from_box.depth / 2)
+
+    t_pt = to_box.frame.point.copy()
+    t_pt -= cg.Point(to_box.width / 2, to_box.height / 2, to_box.depth / 2)
+
+    Tr = cg.Translation.from_vector(t_pt - f_pt)
+
+    R = cg.Rotation.from_frame_to_frame(
+        frame_from=to_box.frame, frame_to=from_box.frame
+    )
+
+    xfactor = from_box.width / to_box.width
+    yfactor =from_box.height / to_box.height
+    zfactor = from_box.depth / to_box.depth
+    S = cg.Scale.from_factors([xfactor, yfactor, zfactor])
+
+    return Tr * R * S

bdm_voxel_builder/helpers/pointcloud.py

@@ -1,3 +1,4 @@
+from collections.abc import Sequence
 import os
 
 import compas.geometry as cg
@@ -36,6 +37,25 @@ def pointcloud_from_ndarray(arr: npt.NDArray, return_values=False):
     return pointcloud
 
 
+def pointcloud_to_grid_array(
+    pointcloud: cg.Pointcloud, grid_size: tuple[int, int, int]
+):
+    """Convert a pointcloud to a grid."""
+    if not isinstance(grid_size, Sequence):
+        grid_size = (grid_size, grid_size, grid_size)
+
+    grid_array = np.zeros(grid_size)
+
+    pts = np.array(pointcloud).round().astype(dtype=int)
+
+    if pts.min() < 0:
+        raise ValueError("Pointcloud contains negative values, needs to be transformed to index grid.") # noqa: E501
+
+    for i, j, k in pts:
+        grid_array[i, j, k] = 1
+    return grid_array
+
+
 def save_pointcloud(
     pointcloud: cg.Pointcloud, values: list[float] = None, note: str = None
 ):

pyproject.toml

@@ -13,6 +13,10 @@ requires-python = ">= 3.11"
 requires = ["setuptools>=61.0"]
 build-backend = "setuptools.build_meta"
 
+[tool.pytest.ini_options]
+addopts = "-vv"
+
+
 [tool.ruff.lint]
 select = [
     # pycodestyle

tests/conftest.py

@@ -0,0 +1,39 @@
+import numpy as np
+import numpy.typing as npt
+import pytest
+
+
+@pytest.fixture
+def random_generator():
+    return np.random.default_rng(1000)
+
+
+@pytest.fixture
+def random_int_array(random_generator):
+    shape = (10, 10, 10)
+    return random_generator.integers(0, 2, size=shape)
+
+
+def _random_pts(n_pts: int, random_generator: np.random.Generator):
+    shape = (n_pts, 3)
+    return 100 * random_generator.random(shape) - 100
+
+@pytest.fixture
+def random_pts():
+    return _random_pts
+
+
+def _activate_random_indices(
+    array: npt.NDArray, random_generator: np.random.Generator
+) -> tuple[npt.NDArray, int]:
+    """Activate random voxels in the grid, return number activated."""
+    random_array = random_generator.integers(0, high=2, size=array.shape).astype(
+        array.dtype
+    )
+
+    return random_array, len(np.flatnonzero(random_array))
+
+
+@pytest.fixture
+def activate_random_indices():
+    return _activate_random_indices

tests/grid/test_base.py

@@ -1,16 +1,18 @@
 import compas.geometry as cg
 import numpy as np
+import numpy.typing as npt
 import pyopenvdb as vdb
 import pytest
 
 from bdm_voxel_builder import get
 from bdm_voxel_builder.grid import Grid
 
-
-@pytest.fixture
-def random_int_array():
-    shape = (10, 10, 10)
-    return np.random.default_rng().integers(0, 2, size=shape)
+def test__activate_random_indices(activate_random_indices ,random_generator):
+    array = np.zeros(shape=(10, 10, 10), dtype=np.float64)
+    array, activated = activate_random_indices(array, random_generator)
+    assert activated == 499
+    assert array[5, 2, 2] == 1
+    assert array[4, 4, 1] == 0
 
 
 def test_to_vdb_grid(random_int_array):
@@ -32,7 +34,7 @@ def test_save_vdb(random_int_array):
 def test_get_local_bbox():
     grid = Grid(grid_size=(10, 10, 10))
     bbox = grid.get_local_bbox()
-    assert bbox.diagonal == ((0, 0, 0), (10, 10, 10))
+    assert bbox.diagonal == cg.Line((0, 0, 0), (10, 10, 10))
 
 
 def test_get_local_bbox_with_non_square_grid_size():
@@ -67,3 +69,96 @@ def test_from_vdb_w_grid():
 
     assert grid.name == "ls_sphere"
     assert grid.array.shape == (124, 124, 124)
+
+
+def test_set_value_at_index():
+    grid = Grid(grid_size=(10, 10, 10))
+    grid.set_value_at_index((5, 5, 5), 1)
+    assert grid.get_value_at_index((5, 5, 5)) == 1
+
+
+def test_get_active_voxels():
+    grid = Grid(grid_size=(10, 10, 10))
+    grid.set_value_at_index((5, 5, 5), 1)
+    active_voxels = grid.get_active_voxels()
+    assert len(active_voxels) == 3
+    assert active_voxels[0][0] == 5
+    assert active_voxels[1][0] == 5
+    assert active_voxels[2][0] == 5
+
+
+def test_get_index_pts(activate_random_indices,random_generator):
+    grid = Grid(grid_size=(10, 10, 10))
+
+    grid.array, activated = activate_random_indices(grid.array, random_generator)
+
+    index_pts = grid.get_index_pts()
+    assert len(index_pts) == activated
+    assert len(index_pts[0]) == 3
+
+
+def test_get_index_pointcloud(activate_random_indices, random_generator):
+    grid = Grid(grid_size=(10, 10, 10))
+
+    grid.array, activated = activate_random_indices(grid.array, random_generator)
+
+    pointcloud = grid.get_index_pointcloud()
+
+    assert len(pointcloud.points) == activated
+
+
+def test_get_world_pts(activate_random_indices, random_generator):
+    grid = Grid(grid_size=(10, 10, 10), xform=cg.Scale.from_factors([2, 2, 2]))
+
+    grid.array, activated = activate_random_indices(grid.array, random_generator)
+
+    world_pts = grid.get_world_pts()
+    assert len(world_pts) == activated
+    assert len(world_pts[0]) == 3
+    assert grid.get_local_bbox().contains_points(world_pts)
+
+
+def test_get_world_pointcloud(activate_random_indices, random_generator):
+    grid = Grid(grid_size=(10, 10, 10), xform=cg.Translation.from_vector([2, 2, 4]))
+
+    grid.array, activated = activate_random_indices(grid.array, random_generator)
+
+    pointcloud = grid.get_world_pointcloud()
+    assert len(pointcloud.points) == activated
+    assert pointcloud.closest_point([0, 0, 0]) == [3, 2, 4]
+
+
+def test_get_world_pointcloud_messy(activate_random_indices, random_generator):
+    grid = Grid(
+        grid_size=(10, 10, 10),
+        xform=cg.Translation.from_vector([2.23, 2.11, 4.13])
+        * cg.Scale.from_factors([1.1, 1.1, 1.1]),
+    )
+
+    grid.array, activated = activate_random_indices(grid.array, random_generator)
+
+    pointcloud = grid.get_world_pointcloud()
+    assert len(pointcloud.points) == activated
+    assert pointcloud.closest_point([0, 0, 0]) == [3.33, 2.11, 4.13]
+
+
+def test_get_merged_array_with():
+    grid1 = Grid(grid_size=(10, 10, 10))
+    grid2 = Grid(grid_size=(10, 10, 10))
+    merged_array = grid1.get_merged_array_with(grid2)
+    assert merged_array.shape == (10, 10, 10)
+
+
+def test_from_pointcloud():
+    pointcloud = cg.Pointcloud([[0, 0, 0], [1, 1, 1], [2, 2, 2]])
+
+    grid = Grid.from_pointcloud(pointcloud, grid_size=3)
+
+    assert grid.array[0, 0, 0] == 1
+    assert grid.array[1, 1, 1] == 1
+    assert grid.array[2, 2, 2] == 1
+
+def test_from_pointcloud_large(random_pts, random_generator):
+    pointcloud = cg.Pointcloud(random_pts(1000, random_generator))
+    grid = Grid.from_pointcloud(pointcloud, grid_size=25)
+    assert grid.get_active_voxels() == 1000