Linear Algebra#
Basics#
Scipp supports basic linear algebra operations in 3 dimensions using vectors of length 3 and several transformations. Variables containing vectors or transformations are created using
scipp.vectors (for an array of vectors)
scipp.vector (for a single vector)
functions in scipp.spatial (for transformations), e.g.,
scipp.spatial.rotations (for and array of rotations)
scipp.spatial.linear_transform (for a single linear transform, i.e. rotation + scaling)
Let us consider an example, creating a 1-D variable containing vectors:
[1]:
import scipp as sc
import scipp.spatial
import numpy as np
import plopp as pp
vecs = sc.vectors(dims=['x'], unit='m', values=np.arange(2 * 3).reshape(2, 3))
vecs
[1]:
- (x: 2)vector3m[0. 1. 2.], [3. 4. 5.]
Values:
array([[0., 1., 2.], [3., 4., 5.]])
While vecs
has only a single dimension x
, the values
property exposes the underlying vector value as a 2-D array of float64
:
[2]:
vecs.values
[2]:
array([[0., 1., 2.],
[3., 4., 5.]])
Access to individual vector components x
, y
, and z
is provided using the special fields
property:
[3]:
vecs.fields.x
[3]:
- (x: 2)float64m0.0, 3.0
Values:
array([0., 3.])
These properties return a view to the underlying data and thus support setting element values as well as in-place modification:
[4]:
vecs.fields.x = 123.0 * sc.units.m
vecs.fields.y += 0.123 * sc.units.m
vecs.values
[4]:
array([[123. , 1.123, 2. ],
[123. , 4.123, 5. ]])
This works in an identical manner for matrices with the properties xx
, xy
, xz
, yx
, yy
, yz
, zx
, zy
, and zz
. The values
property allows for use of, e.g., NumPy functionality. We may, e.g., compute the inverse using numpy.linalg.inv
, which is not supported in Scipp yet:
[5]:
np.random.seed(0) # matrices are not singular for this seed
mats = sc.spatial.linear_transforms(dims=['x'], values=np.random.rand(2, 3, 3))
inv = sc.spatial.linear_transforms(dims=['x'], values=np.linalg.inv(mats.values))
(mats * inv).values.round(decimals=12)
[5]:
array([[[ 1., 0., 0.],
[-0., 1., 0.],
[-0., 0., 1.]],
[[ 1., -0., 0.],
[-0., 1., -0.],
[-0., -0., 1.]]])
Example#
We create some random data and use positions to create a plot three-dimensional scatter plot to illustrate some concepts. Image this as an image sensor with 12x12 pixels:
[6]:
nx = 12
ny = 12
x = sc.linspace(dim='x', start=-0.1, stop=0.1, num=nx, unit='m')
y = sc.linspace(dim='y', start=-0.1, stop=0.1, num=ny, unit='m')
sensor = sc.DataArray(
data=sc.array(dims=['x', 'y'], values=np.random.rand(nx, ny)),
coords={'position': sc.spatial.as_vectors(x, y, 0.0 * sc.units.m)},
)
pp.scatter3d(sensor, pos="position", pixel_size=0.01)
[6]:
We can use the vector-element access properties in combination with slicing to offset some of the pixels:
[7]:
sensor.coords['position']['x', 5:].fields.x += 0.1 * sc.units.m
sensor.coords['position']['y', 5:].fields.y += 0.05 * sc.units.m
pp.scatter3d(sensor, pos="position", pixel_size=0.01)
[7]:
We use sc.spatial.linear_transform
to create a rotation matrix (in this case to rotate by 30 deg around the y
axis) and apply it to the positions:
[8]:
rotation = sc.spatial.linear_transform(
value=[
[0.8660254, 0.0000000, 0.5000000],
[0.0000000, 1.0000000, 0.0000000],
[-0.5000000, 0.0000000, 0.8660254],
]
)
sensor.coords['position']['x', 5:] = rotation * sensor.coords['position']['x', 5:]
pp.scatter3d(sensor, pos="position", pixel_size=0.01)
[8]:
More conveniently, we can create our rotations directly from a rotation vector. Scipp provides the rotations_from_rotvecs
function as part of the scipp.spatial
module. Let’s reverse our original rotation by applying a -30 deg around the y
axis.
[9]:
from scipp.spatial import rotations_from_rotvecs
rotation_back = rotations_from_rotvecs(
rotation_vectors=sc.vector(value=[0, -30.0, 0], unit=sc.units.deg)
)
sensor.coords['position']['x', 5:] = rotation_back * sensor.coords['position']['x', 5:]
pp.scatter3d(sensor, pos="position", pixel_size=0.01)
[9]:
We can compound rotations by multiplying our rotation variable. Lets make rotate the same again in the same axis to give an overall 60 deg around the y
axis, applied to every pixel in our sensor.
[10]:
rotation = rotation * rotation
sensor.coords['position'] = rotation * sensor.coords['position']
pp.scatter3d(sensor, pos="position", pixel_size=0.01)
[10]:
Scipp also provides vector and scalar products of vectors. Lets determine the surface normal vector for our sensor in the new rotation given the knowledge that all points are coplanar.
[11]:
a = sensor.coords['position']['x', 0]['y', 0]
b = sensor.coords['position']['x', -1]['y', 0]
c = sensor.coords['position']['x', 0]['y', -1]
norm = sc.cross(c - a, c - b)
norm /= sc.norm(norm)
norm
[11]:
- ()vector3𝟙[ 0.86602541 -0. 0.5 ]
Values:
array([ 0.86602541, -0. , 0.5 ])
Now lets verify that the normal vector we just calculated corresponds to a rotation of our original normal vector (0, 0, 1) by 60 degrees around y
[12]:
sc.isclose(rotation * sc.vector(value=[0, 0, 1]), norm)
[12]:
- ()boolTrue
Values:
array(True)