Extracting footprints of primitive cubes

I have stages consisting of several primitive cubes. The cubes may be rotated about the up axis, translated, and scaled. The cubes’ up axis scales and translations are the same (i.e., they are the same height and lie on the same plane). I would like a function that extracts the corners of the cubes’ footprints, and a corresponding function that takes a rectangle representing a footprint and reproduces the prim cube.

I have been experimenting with variations of the following, but I have not gotten the translation and/or scaling to work as expected:

def get_transformation_matrix_and_scale(stage, rectangle, group_idx):
    # 1. Determine Translation
    centroid = rectangle.centroid
    translation = [centroid.x, centroid.y, 0]  # Assuming z=0 for 2D space
    
    # 2. Determine Rotation
    # Assuming the first two points form one of the rectangle's edges
    coords = list(rectangle.exterior.coords)
    dx = coords[1][0] - coords[0][0]
    dy = coords[1][1] - coords[0][1]
    angle = np.arctan2(dy, dx)
    
    # 3. Determine Scale
    width = sqrt(dx**2 + dy**2)
    height = sqrt((coords[2][0] - coords[1][0])**2 + (coords[2][1] - coords[1][1])**2)
    scale = [width, height, 1]  # Assuming uniform scaling in z=1
    
    # 4. Compose the Transformation Matrix
    # Rotation matrix around z-axis
    R = np.array([
        [cos(angle), -sin(angle), 0, 0],
        [sin(angle), cos(angle), 0, 0],
        [0, 0, 1, 0],
        [0, 0, 0, 1]
    ])
    print(R)
    
    # Scale matrix
    S = np.diag(scale + [1])
    print(S)
    
    # Translation matrix
    T = np.array([
        [1, 0, 0, translation[0]],
        [0, 1, 0, translation[1]],
        [0, 0, 1, translation[2]],
        [0, 0, 0, 1]
    ])
    print(T)
    
    # Compose the transformation matrix: T * R * S
    transformation_matrix = T @ R @ S


    xform = UsdGeom.Xform.Define(stage, f'/Cube_{group_idx}_grp')
    cube = UsdGeom.Cube.Define(stage, xform.GetPath().AppendChild(f'Cube{group_idx}'))
    transform_op = cube.AddTransformOp(UsdGeom.XformOp.PrecisionDouble)
    transform_op.Set(Gf.Matrix4d(transformation_matrix))
    cube.AddScaleOp(UsdGeom.XformOp.PrecisionDouble).Set(value=Gf.Vec3d(scale))

Please see Linear Algebra in UsdGeom, which states that matrices are row-major with vector per-multiplication, such that translation is specified in the fourth row, not the fourth column. Numpy is also row-major, so the matrix you are specifying is a projection, not a translation.

You could avoid having to keep track of such details by using individual xformOps for translation, rotation, and scale, which would also be a sparser and more easily animatable representation of the transform! :slight_smile:

Sorry for being scatershot here - just noticed your matrix composition, which is also inverted from what you probably want. Because USD pre-multiplies, you likely want S * R * T, not T * R * S

And why are you adding a separate scale op when you have already composed the scale into your composite matrix?

Thank you @spiff ! The order and matrix composition confusion may be my problems.

You could avoid having to keep track of such details by using individual xformOps for translation, rotation, and scale, which would also be a sparser and more easily animatable representation of the transform! :slight_smile:

That sounds good, but I am trying to reproduce the output of a different process that I don’t control. It uses an xformOp:scale and a xformOp:transform for reasons I don’t understand. I think the xformOp:transform represents the rotation and translation only; I did not mean to compose the scale into my transformation matrix.

Thanks again @spiff . I think your answer has gotten me closer, but I don’t quite understand the scale and translation.

I create a stage with 2 cubes as follows. This is meant to represent the data produced by a different, black-box process.

xform = UsdGeom.Xform.Define(stage, f'/Wall_1_grp')
cube = UsdGeom.Cube.Define(stage, xform.GetPath().AppendChild(f'Wall1'))
R = np.eye(4, 4)
s = np.ones(3)
T = np.eye(4, 4)
T[3][0] = 2
T[3][2] = 2
M = Gf.Matrix4d(R @ T)
M.SetRotateOnly(Gf.Rotation(Gf.Vec3d(0, 1, 0), 33))

transform_op = cube.AddTransformOp(UsdGeom.XformOp.PrecisionDouble)
transform_op.Set(M)
cube.AddScaleOp(UsdGeom.XformOp.PrecisionDouble).Set(value=Gf.Vec3d(s[0], s[1], s[2]))

xform2 = UsdGeom.Xform.Define(stage, f'/Wall_2_grp')
cube2 = UsdGeom.Cube.Define(stage, xform2.GetPath().AppendChild(f'Wall2'))
R = np.eye(4, 4)
s = np.ones(3) * 0.5
T = np.eye(4, 4)
T[3][0] = -2
M = Gf.Matrix4d(R @ T)
transform_op2 = cube2.AddTransformOp(UsdGeom.XformOp.PrecisionDouble)
transform_op2.Set(M)
cube2.AddScaleOp(UsdGeom.XformOp.PrecisionDouble).Set(value=Gf.Vec3d(s[0], s[1], s[2]))
stage.GetRootLayer().Save()

Figure 1 visualizes the stage.

I then attempt to parse the .usda and visualize the cubes’ footprints as follows:

import re
import math
import numpy as np
from pxr import Usd, Sdf, Gf, UsdUtils, UsdGeom
from shapely import Polygon
from shapely.affinity import translate, scale, rotate


def is_wall_name(name: str) -> bool:
    return bool(re.compile(r".+Wall\d+$").match(name))
    

stage = Usd.Stage.Open('cubes.usdz')
walls = [prim for prim in stage.Traverse() if is_wall_name(str(prim.GetPath()))]
polygons = []

for wall in walls:
    transform = wall.GetAttribute('xformOp:transform').Get()
    s = np.array(wall.GetAttribute('xformOp:scale').Get())
    R = transform.ExtractRotationMatrix()
    t = transform.ExtractTranslation()
    euler_z_angle = math.atan2(transform[0][1], transform[1][1])
    euler_y_angle = math.atan2(transform[2][0], transform[2][2])
    z_rotation = -(euler_z_angle - euler_y_angle)

    # the side length should be 1...
    unit_square = Polygon([
        [-1, -1],
        [-1, 1],
        [1, 1],
        [1, -1]
    ])

    transformed_polygon = unit_square
    transformed_polygon = scale(transformed_polygon, xfact=s[0], yfact=s[2])
    transformed_polygon = rotate(transformed_polygon, angle=math.degrees(z_rotation))
    transformed_polygon = translate(transformed_polygon, xoff=t[0], yoff=-t[2])
    
    polygons.append(transformed_polygon)


import matplotlib.pyplot as plt
from pylab import rcParams

rcParams['figure.figsize'] = 5, 5
    
fig, ax = plt.subplots()

for i, transformed_polygon in enumerate(polygons):
    x, y = transformed_polygon.exterior.xy
    ax.plot(x, y)
    ax.fill(x, y, alpha=0.4)


ax.set_aspect('equal')

plt.grid()
plt.xlim(-5, 5)
plt.ylim(-5, 5)
plt.show()

This doesn’t look quite right; see figure 2.

If I change the unit_square to have a side length of 2, it looks right (figure 3), but I suspect that is only coincidentally correct.

Am I misunderstanding the prim Cube, or is there an order-of-operations error?

Hi @gavin, I don’t have time to go through your example right now, but I can confirm to you that UsdGeomCube does have a fallback edge/side length of 2. You can change that by authoring, in your scene, both the size and extent properties on each cube, or you can adjust with an xformOp:scale

1 Like

I see, thank you @spiff !