fix: scale VTK mm coordinates to meters when exporting to OpenUSD (#45)

* fix: scale VTK mm coordinates to meters when exporting to OpenUSD

VTK/medical imaging uses millimeters; OpenUSD stages declaring
metersPerUnit=1.0 require meter-scale geometry. Previously ras_points_to_usd()
performed only a RAS→Y-up axis swap, writing raw mm values into USD while the
stage metadata claimed they were meters — making a 100 mm structure appear as
100 meters in viewers.

Changes:
- ras_to_usd() and ras_points_to_usd(): apply * 0.001 (mm → m) during the
  axis-swap so all point coordinates in USD are in meters.
- ras_normals_to_usd(): decoupled from ras_points_to_usd(); performs axis swap
  only with no unit scaling (normals are unit direction vectors).
- usd_tools.py merge_usd_files / merge_usd_files_flattened: fix inconsistent
  metersPerUnit from 0.01 (centimeters) to 1.0 (meters) to match converter output.
- save_usd_file_arrangement: update grid spacing from 400.0 mm to 0.4 m.
- Add TestUnitScaling tests verifying point scaling, normal length preservation,
  and stage metersPerUnit metadata.

Breaking change: all generated USD files will have coordinates 1000× smaller
than before (meter scale instead of millimeter scale).

https://claude.ai/code/session_01L8zowmsyVqXNkZaQivxgc8

* FIX:Fix USD meter scaling and default geometry samples

Author default values for time-sampled mesh points, extents, and normals so
Omniverse/default-time readers can load single-frame and mixed static/animated
USD content correctly.

Allocate float arrays during RAS-to-USD coordinate conversion to avoid
truncating meter-scaled integer inputs, and update USD merge documentation to
match metersPerUnit=1.0.

* DOC: Minor doc update

* COMP: Eliminate copying large array that isn't needed.

---------

Co-authored-by: Claude <noreply@anthropic.com>

Author: aylward

docs/API_MAP.md

@@ -373,19 +373,19 @@ _Re-run `py utils/generate_api_map.py` whenever public APIs change._
 - **class UsdMeshConverter** (line 25): Converts MeshData to UsdGeomMesh with full feature support.
   - `def __init__(self, stage, settings, material_mgr)` (line 36): Initialize mesh converter.
   - `def create_mesh(self, mesh_data, mesh_path, time_code=None, bind_material=True)` (line 53): Create a UsdGeomMesh from MeshData.
-  - `def create_time_varying_mesh(self, mesh_data_sequence, mesh_path, time_codes, bind_material=True)` (line 282): Create a mesh with time-varying attributes.
+  - `def create_time_varying_mesh(self, mesh_data_sequence, mesh_path, time_codes, bind_material=True)` (line 288): Create a mesh with time-varying attributes.
 
 ## src/physiomotion4d/vtk_to_usd/usd_utils.py
 
 - `def ras_to_usd(point)` (line 18): Convert RAS (Right-Anterior-Superior) coordinates to USD's right-handed Y-up system.
-- `def ras_points_to_usd(points)` (line 45): Convert array of RAS points to USD coordinates.
-- `def ras_normals_to_usd(normals)` (line 67): Convert array of RAS normals to USD coordinates.
-- `def numpy_to_vt_array(array, data_type)` (line 81): Convert numpy array to appropriate VtArray type.
-- `def get_sdf_value_type(data_type, num_components)` (line 153): Get appropriate SDF value type for primvar creation.
-- `def sanitize_primvar_name(name)` (line 200): Sanitize a name to be USD-compliant.
-- `def create_primvar(geom, array, array_name_prefix='', time_code=None)` (line 235): Create a USD primvar from a GenericArray.
-- `def triangulate_face(face_counts, face_indices)` (line 349): Triangulate polygonal faces.
-- `def compute_mesh_extent(points)` (line 389): Compute bounding box extent for a mesh.
+- `def ras_points_to_usd(points)` (line 53): Convert array of RAS points (mm) to USD coordinates (m).
+- `def ras_normals_to_usd(normals)` (line 76): Convert array of RAS normals to USD Y-up coordinates.
+- `def numpy_to_vt_array(array, data_type)` (line 99): Convert numpy array to appropriate VtArray type.
+- `def get_sdf_value_type(data_type, num_components)` (line 171): Get appropriate SDF value type for primvar creation.
+- `def sanitize_primvar_name(name)` (line 218): Sanitize a name to be USD-compliant.
+- `def create_primvar(geom, array, array_name_prefix='', time_code=None)` (line 253): Create a USD primvar from a GenericArray.
+- `def triangulate_face(face_counts, face_indices)` (line 367): Triangulate polygonal faces.
+- `def compute_mesh_extent(points)` (line 407): Compute bounding box extent for a mesh.
 
 ## src/physiomotion4d/vtk_to_usd/vtk_reader.py
 
@@ -741,6 +741,10 @@ _Re-run `py utils/generate_api_map.py` whenever public APIs change._
   - `def test_time_series_conversion(self, test_directories, kcl_average_surface)` (line 517): Test converting multiple VTK files as a time series.
 - **class TestIntegration** (line 557): Integration tests combining multiple features.
   - `def test_end_to_end_conversion(self, test_directories, kcl_average_surface)` (line 560): Test complete conversion workflow with all features.
+- **class TestUnitScaling** (line 603): Verify that VTK mm coordinates are converted to USD meter coordinates.
+  - `def test_mm_to_m_point_scaling(self, tmp_path)` (line 606): Points written to USD must be 0.001× their original mm values.
+  - `def test_normals_remain_unit_length(self, tmp_path)` (line 638): Normal vectors must not be scaled — they should remain unit length.
+  - `def test_stage_meters_per_unit(self, tmp_path)` (line 664): Stage metersPerUnit metadata must be 1.0 (coordinates stored in meters).
 
 ## utils/claude_github_reviews.py
 

src/physiomotion4d/usd_tools.py

@@ -170,8 +170,8 @@ def save_usd_file_arrangement(
         n_rows = int(np.floor(np.sqrt(n_objects)))
         n_cols = int(np.ceil(n_objects / n_rows))
         self.log_info("Grid layout: %d rows x %d cols", n_rows, n_cols)
-        x_spacing = 400.0
-        y_spacing = 400.0
+        x_spacing = 0.4
+        y_spacing = 0.4
         x_offset = -x_spacing * (n_cols - 1) / 2
         y_offset = -y_spacing * (n_rows - 1) / 2
 
@@ -273,8 +273,8 @@ def merge_usd_files(
                 coordinate systems
 
         Note:
-            The merged file uses meters as the base unit (0.01 scale factor)
-            and Y-up axis orientation, which are standard for Omniverse.
+            The merged file stores coordinates in meters (metersPerUnit=1.0)
+            with upAxis="Y", which are standard for Omniverse.
             Time-varying data (animations) are preserved across all time samples.
 
         Example:
@@ -285,7 +285,7 @@ def merge_usd_files(
         """
         # Create new stage with meters as units (standard USD configuration)
         stage = Usd.Stage.CreateNew(output_filename)
-        stage.SetMetadata("metersPerUnit", 0.01)
+        stage.SetMetadata("metersPerUnit", 1.0)
         stage.SetMetadata("upAxis", "Y")
 
         # Define root prim for organization
@@ -478,7 +478,7 @@ def merge_usd_files_flattened(
         temp_stage = Usd.Stage.CreateInMemory()
 
         # Set standard metadata (meters and Y-up for Omniverse)
-        temp_stage.SetMetadata("metersPerUnit", 0.01)
+        temp_stage.SetMetadata("metersPerUnit", 1.0)
         temp_stage.SetMetadata("upAxis", "Y")
 
         # Define root prim for organization

src/physiomotion4d/vtk_to_usd/usd_mesh_converter.py

@@ -97,6 +97,8 @@ def create_mesh(
         # Set points (time-varying if time_code provided)
         points_attr = mesh.CreatePointsAttr()
         if time_code is not None:
+            if points_attr.Get() is None:
+                points_attr.Set(usd_points)
             points_attr.Set(usd_points, time_code)
         else:
             points_attr.Set(usd_points)
@@ -105,6 +107,8 @@ def create_mesh(
         extent = compute_mesh_extent(usd_points)
         extent_attr = mesh.CreateExtentAttr()
         if time_code is not None:
+            if extent_attr.Get() is None:
+                extent_attr.Set(extent)
             extent_attr.Set(extent, time_code)
         else:
             extent_attr.Set(extent)
@@ -120,6 +124,8 @@ def create_mesh(
             normals_attr = mesh.CreateNormalsAttr()
             normals_attr.SetMetadata("interpolation", UsdGeom.Tokens.vertex)
             if time_code is not None:
+                if normals_attr.Get() is None:
+                    normals_attr.Set(usd_normals)
                 normals_attr.Set(usd_normals, time_code)
             else:
                 normals_attr.Set(usd_normals)

src/physiomotion4d/vtk_to_usd/usd_utils.py

@@ -28,54 +28,72 @@ def ras_to_usd(point: NDArray | tuple | list) -> Gf.Vec3f:
     - Y: up
     - Z: back (toward camera)
 
-    Conversion: USD(x, y, z) = RAS(x, z, -y)
+    Conversion: USD(x, y, z) = RAS(x, z, -y) * 0.001  (mm → m)
 
     Args:
-        point: Point in RAS coordinates [x, y, z]
+        point: Point in RAS coordinates [x, y, z] in millimeters
 
     Returns:
-        Gf.Vec3f: Point in USD coordinates
+        Gf.Vec3f: Point in USD coordinates in meters
     """
     if isinstance(point, (tuple, list)):
-        return Gf.Vec3f(float(point[0]), float(point[2]), float(-point[1]))
+        return Gf.Vec3f(
+            float(point[0]) * 0.001,
+            float(point[2]) * 0.001,
+            float(-point[1]) * 0.001,
+        )
     else:
-        return Gf.Vec3f(float(point[0]), float(point[2]), float(-point[1]))
+        return Gf.Vec3f(
+            float(point[0]) * 0.001,
+            float(point[2]) * 0.001,
+            float(-point[1]) * 0.001,
+        )
 
 
 def ras_points_to_usd(points: NDArray) -> Vt.Vec3fArray:
-    """Convert array of RAS points to USD coordinates.
+    """Convert array of RAS points (mm) to USD coordinates (m).
+
+    Applies axis swap RAS → Y-up and scales millimeters to meters (* 0.001).
 
     Args:
-        points: Array of points with shape (N, 3)
+        points: Array of points with shape (N, 3) in millimeters
 
     Returns:
-        Vt.Vec3fArray: Points in USD coordinates
+        Vt.Vec3fArray: Points in USD Y-up coordinates in meters
     """
     if points.shape[1] != 3:
         raise ValueError(f"Points must have shape (N, 3), got {points.shape}")
 
-    # Vectorized conversion: USD(x, y, z) = RAS(x, z, -y)
-    usd_points = np.empty_like(points)
-    usd_points[:, 0] = points[:, 0]  # X stays the same
-    usd_points[:, 1] = points[:, 2]  # Y = Z
-    usd_points[:, 2] = -points[:, 1]  # Z = -Y
+    # Vectorized: USD(x, y, z) = RAS(x, z, -y) * 0.001  (mm → m)
+    usd_points = np.empty(points.shape, dtype=np.float32)
+    usd_points[:, 0] = points[:, 0] * 0.001
+    usd_points[:, 1] = points[:, 2] * 0.001
+    usd_points[:, 2] = -points[:, 1] * 0.001
 
-    # Convert to USD Vec3fArray
-    return Vt.Vec3fArray.FromNumpy(usd_points.astype(np.float32))
+    return Vt.Vec3fArray.FromNumpy(usd_points)
 
 
 def ras_normals_to_usd(normals: NDArray) -> Vt.Vec3fArray:
-    """Convert array of RAS normals to USD coordinates.
+    """Convert array of RAS normals to USD Y-up coordinates.
 
-    Same transformation as points since normals are vectors.
+    Applies only the axis swap — normals are unit direction vectors and must
+    not be scaled by the mm→m factor.
 
     Args:
         normals: Array of normals with shape (N, 3)
 
     Returns:
-        Vt.Vec3fArray: Normals in USD coordinates
+        Vt.Vec3fArray: Normals in USD Y-up coordinates (unit length preserved)
     """
-    return ras_points_to_usd(normals)
+    if normals.shape[1] != 3:
+        raise ValueError(f"Normals must have shape (N, 3), got {normals.shape}")
+
+    usd_normals = np.empty(normals.shape, dtype=np.float32)
+    usd_normals[:, 0] = normals[:, 0]
+    usd_normals[:, 1] = normals[:, 2]
+    usd_normals[:, 2] = -normals[:, 1]
+
+    return Vt.Vec3fArray.FromNumpy(usd_normals)
 
 
 def numpy_to_vt_array(array: NDArray, data_type: DataType) -> Any:

tests/test_vtk_to_usd_library.py

@@ -598,3 +598,77 @@ def test_end_to_end_conversion(
         print(f"  Size: {output_usd.stat().st_size / 1024:.1f} KB")
         print(f"  Points: {len(points):,}")
         print(f"  Primvars: {len(primvars)}")
+
+
+class TestUnitScaling:
+    """Verify that VTK mm coordinates are converted to USD meter coordinates."""
+
+    def test_mm_to_m_point_scaling(self, tmp_path: Path) -> None:
+        """Points written to USD must be 0.001× their original mm values."""
+        # Sphere with radius=100 mm — vertices should be near ±100 in VTK.
+        mesh = pv.Sphere(radius=100.0)
+        output_usd = tmp_path / "sphere.usd"
+
+        stage = ConvertVTKToUSD(
+            data_basename="Sphere",
+            input_polydata=[mesh],
+        ).convert(str(output_usd))
+
+        mesh_prim = stage.GetPrimAtPath("/World/Sphere/Mesh")
+        assert mesh_prim.IsValid(), "Mesh prim not found at expected path"
+
+        usd_mesh = UsdGeom.Mesh(mesh_prim)
+        usd_points = usd_mesh.GetPointsAttr().Get()
+        assert usd_points is not None and len(usd_points) > 0
+
+        coords = np.array(usd_points)
+        max_coord = float(np.abs(coords).max())
+
+        # In meters a 100 mm sphere has vertices ≤ 0.1 m (plus floating-point headroom).
+        assert max_coord < 0.15, (
+            f"Max coordinate {max_coord:.4f} is not in meters. "
+            "Expected < 0.15 m for a 100 mm radius sphere; "
+            "got a value that looks like millimeters."
+        )
+        # Sanity-check it's not collapsed to near zero (e.g., double-scaling).
+        assert max_coord > 0.05, (
+            f"Max coordinate {max_coord:.6f} is unexpectedly small."
+        )
+
+    def test_normals_remain_unit_length(self, tmp_path: Path) -> None:
+        """Normal vectors must not be scaled — they should remain unit length."""
+        mesh = pv.Sphere(radius=100.0)
+        mesh.compute_normals(inplace=True)
+        output_usd = tmp_path / "sphere_normals.usd"
+
+        stage = ConvertVTKToUSD(
+            data_basename="Sphere",
+            input_polydata=[mesh],
+        ).convert(str(output_usd))
+
+        mesh_prim = stage.GetPrimAtPath("/World/Sphere/Mesh")
+        usd_mesh = UsdGeom.Mesh(mesh_prim)
+        normals_attr = usd_mesh.GetNormalsAttr()
+
+        if normals_attr is None or normals_attr.Get() is None:
+            pytest.skip("No normals on this mesh")
+
+        normals = np.array(normals_attr.Get())
+        norms = np.linalg.norm(normals, axis=1)
+        # Every normal should be ≈ 1.0 (unit vector), not 0.001.
+        assert np.allclose(norms, 1.0, atol=1e-3), (
+            f"Normals are not unit length after conversion. "
+            f"Mean norm: {norms.mean():.6f}, expected ≈ 1.0"
+        )
+
+    def test_stage_meters_per_unit(self, tmp_path: Path) -> None:
+        """Stage metersPerUnit metadata must be 1.0 (coordinates stored in meters)."""
+        mesh = pv.Sphere(radius=100.0)
+        output_usd = tmp_path / "sphere_meta.usd"
+
+        stage = ConvertVTKToUSD(
+            data_basename="Sphere",
+            input_polydata=[mesh],
+        ).convert(str(output_usd))
+
+        assert UsdGeom.GetStageMetersPerUnit(stage) == 1.0