Menu
I have a voxel (np.array) with size 3x3x3, filled with some values, this setup is essential for me. I want to have rotation-invariant representation of it. For this case, I decided to try PCA representation which is believed to be invariant to orthogonal transformations. another
For simplicity, I took some axes swap, but in case I'm mistaken there can be np.rot90.
I have interpereted my 3d voxels as a set of weighted 3d cube point vectors which I incorrectly called "basis", total 27 (so that is some set of 3d point in space, represented by the vectors, obtained from cube points, scaled by voxel values).
import numpy as np
voxel1 = np.random.normal(size=(3,3,3))
voxel2 = np.transpose(voxel1, (1,0,2)) #np.rot90(voxel1) #
basis = []
for i in range(3):
for j in range(3):
for k in range(3):
basis.append([i+1, j+1, k+1]) # avoid 0
basis = np.array(basis)
voxel1 = voxel1.reshape((27,1))
voxel2 = voxel2.reshape((27,1))
voxel1 = voxel1*basis # weighted basis vectors
voxel2 = voxel2*basis
print(voxel1.shape)
(27, 3)
Then I did PCA to those 27 3-dimensional vectors:
def pca(x):
center = np.mean(x, 0)
x = x - center
cov = np.cov(x.T) / x.shape[0]
e_values, e_vectors = np.linalg.eig(cov)
order = np.argsort(e_values)
v = e_vectors[:, order].transpose()
return x.dot(v)
vp1 = pca(voxel1)
vp2 = pca(voxel2)
But the results in vp1 and vp2 are different. Perhaps, I have a mistake (though I beleive this is the right formula), and the proper code must be
x.dot(v.T)
But in this case the results are very strange. The upper and bottom blocks of the transofrmed data are the same up to the sign:
>>> np.abs(np.abs(vp1)-np.abs(vp2)) > 0.01
array([[False, False, False],
[False, False, False],
[False, False, False],
[ True, True, True],
[ True, True, True],
[ True, True, True],
[ True, True, True],
[ True, True, True],
[ True, False, True],
[ True, True, True],
[ True, True, True],
[ True, True, True],
[False, False, False],
[False, False, False],
[False, False, False],
[ True, True, True],
[ True, True, True],
[ True, True, True],
[ True, True, True],
[ True, True, True],
[ True, False, True],
[ True, True, True],
[ True, True, True],
[ True, True, True],
[False, False, False],
[False, False, False],
[False, False, False]])
What I'm doing wrong?
What I want to do is to find some invariant representation of my weighted voxel, something like positioning according to the axes of inertia or principal axes. I would really appreciate if someone helps me.
UPD: Found the question similar to mine, but code is unavailable
EDIT2: Found the code InertiaRotate and managed to monkey-do the following:
import numpy as np
# https://github.com/smparker/orient-molecule/blob/master/orient.py
voxel1 = np.random.normal(size=(3,3,3))
voxel2 = np.transpose(voxel1, (1,0,2))
voxel1 = voxel1.reshape((27,))
voxel2 = voxel2.reshape((27,))
basis = []
for i in range(3):
for j in range(3):
for k in range(3):
basis.append([i+1, j+1, k+1]) # avoid 0
basis = np.array(basis)
basis = basis - np.mean(basis, axis=0)
def rotate_func(data, mass):
#mass = [ masses[n.lower()] for n in geom.names ]
inertial_tensor = -np.einsum("ax,a,ay->xy", data, mass, data)
# negate sign to reverse the sorting of the tensor
eig, axes = np.linalg.eigh(-inertial_tensor)
axes = axes.T
# adjust sign of axes so third moment moment is positive new in X, and Y axes
testcoords = np.dot(data, axes.T) # a little wasteful, but fine for now
thirdmoment = np.einsum("ax,a->x", testcoords**3, mass)
for i in range(2):
if thirdmoment[i] < 1.0e-6:
axes[i,:] *= -1.0
# rotation matrix must have determinant of 1
if np.linalg.det(axes) < 0.0:
axes[2,:] *= -1.0
return axes
axes1 = rotate_func(basis, voxel1)
v1 = np.dot(basis, axes1.T)
axes2 = rotate_func(basis, voxel2)
v2 = np.dot(basis, axes2.T)
print(v1)
print(v2)
It seems to use basis (coordinates) and mass separately. The results are quite similar to my problem above: some parts of the transformed data match up to the sign, I believe those are some cube sides
print(np.abs(np.abs(v1)-np.abs(v2)) > 0.01)
[[False False False]
[False False False]
[False False False]
[ True True True]
[ True True True]
[ True True True]
[ True True True]
[False False False]
[ True True True]
[ True True True]
[ True True True]
[ True True True]
[False False False]
[False False False]
[False False False]
[ True True True]
[ True True True]
[ True True True]
[ True True True]
[False False False]
[ True True True]
[ True True True]
[ True True True]
[ True True True]
[False False False]
[False False False]
[False False False]]
Looking for some explanation. This code is designed for molecules, and must work...
UPD: Tried to choose 3 vectors as a new basis from those 24 - the one with biggest norm, the one with the smallest and their cross product. Combined them into the matrix V, then used the formula V^(-1)*X to transform coordinates, and got the same problem - the resulting sets of vectors are not equal for rotated voxels.
UPD2: I agree with meTchaikovsky that my idea of multiplying voxel vectors by weights and thus creating some non-cubic point cloud was incorrect. Probably, we indeed need to take the solution for rotated "basis"(yes, this is not a basis, but rather a way to determine point cloud) which will work later when "basis" is the same, but the weights are rotated according to the 3D rotation.
Based on the answer and the reference provided by meTchaikovsky, and finding other answers we together with my friend came to conclusion that rotate_func from molecular package mentioned above tries to invent some convention for computing the signs of the components. Their solution tries to use 3rd moment for the first 2 axes and determinant for the last axis (?). We tried a bit another approach and succeeded to have half of the representations matching:
# -*- coding: utf-8 -*-
"""
Created on Fri Oct 16 11:40:30 2020
@author: Dima
"""
import numpy as np
from numpy.random import randn
from numpy import linalg as la
from scipy.spatial.transform import Rotation as R
np.random.seed(10)
rotate_mat = lambda theta: np.array([[np.cos(theta),-np.sin(theta),0.],[np.sin(theta),np.cos(theta),0.],[0.,0.,1.]])
def pca(feat, x):
# pca with attemt to create convention on sign changes
x_c =x- np.mean(x,axis=0)
x_f= feat*x
x_f-= np.mean(x_f, axis=0)
cov = np.cov(x_f.T)
e_values, e_vectors = np.linalg.eig(cov)
order = np.argsort(e_values)[::-1]
#print(order)
v = e_vectors[:,order]
v= v/np.sign(v[0,:])
if(la.det(v)<0):
v= -v
return x_c @ v
def standardize(x):
# take vector with biggest norm, with smallest and thir cross product as basis
x -= np.mean(x,axis=0)
nrms= la.norm(x, axis=1)
imin= argmin(nrms)
imax= argmax(nrms)
vec1= x[imin, :]
vec2= x[imax, :]
vec3= np.cross(vec1, vec2)
Smat= np.stack([vec1, vec2, vec3], axis=0)
if(la.det(Smat)<0):
Smat= -Smat
return(la.inv(Smat)@x.T)
angles = np.linspace(0.0,90.0,91)
voxel1 = np.random.normal(size=(3,3,3))
res = []
for angle in angles:
voxel2 = voxel1.copy()
voxel1 = voxel1.reshape(27,1)
voxel2 = voxel2.reshape(27,1)
basis1 = np.array([[i+1,j+1,k+1] for k in range(3) for j in range(3) for i in range(3)]).astype(np.double)
basis1 = basis1+1e-4*randn(27,3) # perturbation
basis2 = basis1 @rotate_mat(np.deg2rad(angle))
#voxel1 = voxel1*basis1
#voxel2 = voxel2*basis2
#print(angle,(np.abs(pca(voxel1) - pca(voxel2) )))
#gg= np.abs(standardize(basis1) - standardize(basis2) )
gg= np.abs(pca(voxel1, basis1) - pca(voxel1, basis2) )
ss= np.sum(np.ravel(gg))
bl= np.all(gg<1e-4)
print(angle,ss, bl)
#res.append(np.all(np.abs(pca(voxel1) - pca(voxel2) < 1e-6)))
del basis1, basis2
The results are good up to 58 degree angle (yet we're still experimenting with rotation of x, y axes). After that we have constant difference which indicates some uncounted sign reverse. This is better than the less consistent result of rotate_func:
0.0 0.0 True
1.0 1.1103280567106161e-13 True
2.0 5.150139890290964e-14 True
3.0 8.977126225544196e-14 True
4.0 5.57341699240722e-14 True
5.0 4.205149954378956e-14 True
6.0 3.7435437643664957e-14 True
7.0 1.2943967187158123e-13 True
8.0 5.400185371573149e-14 True
9.0 8.006410204958181e-14 True
10.0 7.777189536904011e-14 True
11.0 5.992073021576436e-14 True
12.0 6.3716122222085e-14 True
13.0 1.0120048110065158e-13 True
14.0 1.4193029076233626e-13 True
15.0 5.32774440341853e-14 True
16.0 4.056702432878251e-14 True
17.0 6.52062429116855e-14 True
18.0 1.3237663595853556e-13 True
19.0 8.950259695710006e-14 True
20.0 1.3795067925438317e-13 True
21.0 7.498727794307339e-14 True
22.0 8.570866862371226e-14 True
23.0 8.961510590826412e-14 True
24.0 1.1839169916779899e-13 True
25.0 1.422193407555868e-13 True
26.0 6.578778015788652e-14 True
27.0 1.0042963537887101e-13 True
28.0 8.438153062569065e-14 True
29.0 1.1299103064863272e-13 True
30.0 8.192453876745831e-14 True
31.0 1.2618492405483406e-13 True
32.0 4.9237819394886296e-14 True
33.0 1.0971028569666842e-13 True
34.0 1.332138304559801e-13 True
35.0 5.280024600049296e-14 True
From the code above, you can see that we tried to use another basis: vector with the biggest norm, vector with the smallest and their cross product. Here we should have only two variants (direction of the cross product) which could be later fixed, but I couldn't manage this alternative solution to work.
I hope that someone can help me finish this and obtain rotation-invariant representation for voxels.
EDIT 3. Thank you very much meTchaikovsky, but the situation is still unclear. My problem initially lies in processing 3d voxels which are (3,3,3) numpy arrays. We reached the conclusion that for finding invariant representation, we just need to fix 3d voxel as weights used for calculating cov matrix, and apply rotations on the centered "basis" (some vectors used for describing point cloud).
Therefore, when we achieved invariance to "basis" rotations, the problem should have been solved: now, when we fix "basis" and use rotated voxel, the result must be invariant. Surprisingly, this is not so. Here I check 24 rotations of the cube with basis2=basis1 (except small perturbation):
import scipy.ndimage
def pca(feat, x):
# pca with attemt to create convention on sign changes
x_c = x - np.mean(x,axis=0)
x_f = feat * x
x_f -= np.mean(x_f,axis=0)
cov = np.cov(x_f.T)
e_values, e_vectors = np.linalg.eig(cov)
order = np.argsort(e_values)[::-1]
v = e_vectors[:,order]
# here is the solution, we switch the sign of the projections
# so that the projection with the largest absolute value along a principal axis is positive
proj = x_c @ v
asign = np.sign(proj)
max_ind = np.argmax(np.abs(proj),axis=0)[None,:]
sign = np.take_along_axis(asign,max_ind,axis=0)
proj = proj * sign
return proj
def rotate_3d(image1, alpha, beta, gamma):
# z
# The rotation angle in degrees.
image2 = scipy.ndimage.rotate(image1, alpha, mode='nearest', axes=(0, 1), reshape=False)
# rotate along y-axis
image3 = scipy.ndimage.rotate(image2, beta, mode='nearest', axes=(0, 2), reshape=False)
# rotate along x-axis
image4 = scipy.ndimage.rotate(image3, gamma, mode='nearest', axes=(1, 2), reshape=False)
return image4
voxel10 = np.random.normal(size=(3,3,3))
angles = [[x,y,z] for x in [-90,0,90] for y in [-90,0,90] for z in [-90,0,90]]
res = []
for angle in angles:
voxel2 = rotate_3d(voxel10, angle[0], angle[1], angle[2])
voxel1 = voxel10.reshape(27,1)
voxel2 = voxel2.reshape(27,1)
basis1 = np.array([[i+1,j+1,k+1] for k in range(3) for j in range(3) for i in range(3)]).astype(np.double)
basis1 += 1e-4*np.random.normal(size=(27, 1)) # perturbation
basis2 = basis1
original_diff = np.sum(np.abs(basis1-basis2))
gg= np.abs(pca(voxel1, basis1) - pca(voxel2, basis2))
ss= np.sum(np.ravel(gg))
bl= np.all(gg<1e-4)
print('difference before pca %.3f,' % original_diff, 'difference after pca %.3f' % ss,bl)
res.append(bl)
del basis1, basis2
print('correct for %.1f percent of time' % (100*(np.sum(res) / len(res))))
difference before pca 0.000, difference after pca 45.738 False
difference before pca 0.000, difference after pca 12.157 False
difference before pca 0.000, difference after pca 26.257 False
difference before pca 0.000, difference after pca 37.128 False
difference before pca 0.000, difference after pca 52.131 False
difference before pca 0.000, difference after pca 45.436 False
difference before pca 0.000, difference after pca 42.226 False
difference before pca 0.000, difference after pca 18.959 False
difference before pca 0.000, difference after pca 38.888 False
difference before pca 0.000, difference after pca 12.157 False
difference before pca 0.000, difference after pca 26.257 False
difference before pca 0.000, difference after pca 50.613 False
difference before pca 0.000, difference after pca 52.132 False
difference before pca 0.000, difference after pca 0.000 True
difference before pca 0.000, difference after pca 52.299 False
Here basis1=basis2 (hence basis difference before pca=0), and you can see 0 for (0,0,0) rotation. But rotated voxels give different result. In case scipy does something wrong, I've checked the approach with numpy.rot90 with the same result:
rot90 = np.rot90
def rotations24(polycube):
# imagine shape is pointing in axis 0 (up)
# 4 rotations about axis 0
yield from rotations4(polycube, 0)
# rotate 180 about axis 1, now shape is pointing down in axis 0
# 4 rotations about axis 0
yield from rotations4(rot90(polycube, 2, axis=1), 0)
# rotate 90 or 270 about axis 1, now shape is pointing in axis 2
# 8 rotations about axis 2
yield from rotations4(rot90(polycube, axis=1), 2)
yield from rotations4(rot90(polycube, -1, axis=1), 2)
# rotate about axis 2, now shape is pointing in axis 1
# 8 rotations about axis 1
yield from rotations4(rot90(polycube, axis=2), 1)
yield from rotations4(rot90(polycube, -1, axis=2), 1)
def rotations4(polycube, axis):
"""List the four rotations of the given cube about the given axis."""
for i in range(4):
yield rot90(polycube, i, axis)
def rot90(m, k=1, axis=2):
"""Rotate an array k*90 degrees in the counter-clockwise direction around the given axis"""
m = np.swapaxes(m, 2, axis)
m = np.rot90(m, k)
m = np.swapaxes(m, 2, axis)
return m
voxel10 = np.random.normal(size=(3,3,3))
gen = rotations24(voxel10)
res = []
for voxel2 in gen:
#voxel2 = rotate_3d(voxel10, angle[0], angle[1], angle[2])
voxel1 = voxel10.reshape(27,1)
voxel2 = voxel2.reshape(27,1)
basis1 = np.array([[i+1,j+1,k+1] for k in range(3) for j in range(3) for i in range(3)]).astype(np.double)
basis1 += 1e-4*np.random.normal(size=(27, 1)) # perturbation
basis2 = basis1
original_diff = np.sum(np.abs(basis1-basis2))
gg= np.abs(pca(voxel1, basis1) - pca(voxel2, basis2))
ss= np.sum(np.ravel(gg))
bl= np.all(gg<1e-4)
print('difference before pca %.3f,' % original_diff, 'difference after pca %.3f' % ss,bl)
res.append(bl)
del basis1, basis2
print('correct for %.1f percent of time' % (100*(np.sum(res) / len(res))))
I tried to investigate this case, and the only perhaps irrelevant thing I found the following:
voxel1 = np.ones((3,3,3))
voxel1[0,0,0] = 0 # if I change 0 to 0.5 it stops working at all
# mirrored around diagonal
voxel2 = np.ones((3,3,3))
voxel2[2,2,2] = 0
for angle in range(1):
voxel1 = voxel1.reshape(27,1)
voxel2 = voxel2.reshape(27,1)
basis1 = np.array([[i+1,j+1,k+1] for k in range(3) for j in range(3) for i in range(3)]).astype(np.double)
basis1 = basis1 + 1e-4 * randn(27,3) # perturbation
basis2 = basis1
# If perturbation is used we have
# difference before pca 0.000, difference after pca 0.000 True
# correct for 100.0 percent of time
# eigenvalues for both voxels
# [1.03417495 0.69231107 0.69235402]
# [0.99995368 0.69231107 0.69235402]
# If no perturbation applied for basis, difference is present
# difference before pca 0.000, difference after pca 55.218 False
# correct for 0.0 percent of time
# eignevalues for both voxels (always have 1.):
# [0.69230769 1.03418803 0.69230769]
# [1. 0.69230769 0.69230769]
Currently don't know how to proceed from there.
EDIT4:
I'm currently thinking that there is some problem with voxel rotations transformed into basis coefficients via voxel.reshape()
Simple experiment with creating array of indices
indices = np.arange(27)
indices3d = indices.reshape((3,3,3))
voxel10 = np.random.normal(size=(3,3,3))
assert voxel10[0,1,2] == voxel10.ravel()[indices3d[0,1,2]]
And then using it for rotations
gen = rotations24(indices3d)
res = []
for ind2 in gen:
basis1 = np.array([[i+1,j+1,k+1] for k in range(3) for j in range(3) for i in range(3)]).astype(np.double)
voxel1 = voxel10.copy().reshape(27,1) #np.array([voxel10[i,j,k] for k in range(3) for j in range(3) for i in range(3)])[...,np.newaxis]
voxel2 = voxel1[ind2.reshape(27,)]
basis1 += 1e-4*np.random.normal(size=(27, 1)) # perturbation
basis2 = basis1[ind2.reshape(27,)]
original_diff = np.sum(np.abs(basis1-basis2))
gg= np.abs(pca(voxel1, basis1) - pca(voxel2, basis2))
ss= np.sum(np.ravel(gg))
bl= np.all(gg<1e-4)
print('difference before pca %.3f,' % original_diff, 'difference after pca %.3f' % ss,bl)
res.append(bl)
del basis1, basis2
print('correct for %.1f percent of time' % (100*(np.sum(res) / len(res))))
Shows that those rotations are not correct, because on my opinion rotated voxel and basis should match:
difference before pca 0.000, difference after pca 0.000 True
difference before pca 48.006, difference after pca 87.459 False
difference before pca 72.004, difference after pca 70.644 False
difference before pca 48.003, difference after pca 71.930 False
difference before pca 72.004, difference after pca 79.409 False
difference before pca 84.005, difference after pca 36.177 False
EDIT 5: Okaaay, so here we go at least for 24 rotations. At first, we had a slight change of logic lurked into our pca function. Here we center x_c (basis) and forget about it, further centering x_f (features*basis) and transforming it with pca. This does not work perhaps because our basis is not centered and multiplication by features further increased the bias. If we center x_c first, and multiply it by features, everything will be Ok. Also, previously we had proj = x_c @ v with v computed from x_f which was totally wrong in this case, as x_f and x_c were centered around different centers.
def pca(feat, x):
# pca with attemt to create convention on sign changes
x_c = x - np.mean(x,axis=0)
x_f = feat * x
x_f -= np.mean(x_f,axis=0)
cov = np.cov(x_f.T)
e_values, e_vectors = np.linalg.eig(cov)
order = np.argsort(e_values)[::-1]
v = e_vectors[:,order]
# here is the solution, we switch the sign of the projections
# so that the projection with the largest absolute value along a principal axis is positive
proj = x_f @ v
return proj
Secondly, as we already found, we need to sort vectors obtained by pca, for example by the first column:
basis2 = basis1
original_diff = np.sum(np.abs(basis1-basis2))
a = pca(voxel1, basis1)
t1 = a[a[:,0].argsort()]
a = pca(voxel2, basis2)
t2 = a[a[:,0].argsort()]
gg= np.abs(t1-t2)
And the last thing we also discovered already, is that simple reshape is wrong for voxel, it must correspond to rotation:
voxel2 = voxel1[ind2.reshape(27,)] #np.take(voxel10, ind2).reshape(27,1).
One more important comment to understand the solution. When we perform PCA on the 3d vectors (point cloud, defined by our basis) with weights assigned (analogously to the inertia of the rigid body), the actual assignment of the weights to the points is sort of external information, which becomes hard-defined for the algorithm. When we rotated basis by applying rotation matrices, we did not change the order of the vectors in the array, hence the order of the mass assignments wasn't changed too. When we start to rotate voxel, we change the order of the masses, so in general PCA algorithm will not work without the same transformation applied to the basis. So, only if we have some array of 3d vectors, transformed by some rotation AND the list of masses re-arranged accordingly, we can detect the rotation of the rigid body using PCA. Otherwise, if we detach masses from points, that would be another body in general.
So how does it work for us then? It works because our points are fully symmetric around the center after centering basis. In this case reassignment of the masses does not change "the body" because vector norms are the same. In this case we can use the same (numerically) basis2=basis1 for testing 24 rotations and rotated voxel2 (rotated point cloud cubes match, just masses migrate). This correspond to the rotation of the point cloud with mass points around the center of the cube. PCA will transform vectors with the same lengths and different masses in the same way according to the body's "inertia" then (after we reached convention on the signs of the components). The only thing left is to sort the pca transformed vectors in the end, because they have different position in the array (because our body was rotated, mass points changed their positions). This makes us lose some information related to the order of the vectors but it looks inevitable.
Here is the code which checks the solution for 24 rotations. If should theoretically work in the general case as well, giving some closer values for more complicated objects rotated inside a bigger voxel:
import numpy as np
from numpy.random import randn
#np.random.seed(20)
def pca(feat, x):
# pca with attemt to create convention on sign changes
x_c = x - np.mean(x,axis=0)
x_f = feat * x_c
cov = np.cov(x_f.T)
e_values, e_vectors = np.linalg.eig(cov)
order = np.argsort(e_values)[::-1]
v = e_vectors[:,order]
# here is the solution, we switch the sign of the projections
# so that the projection with the largest absolute value along a principal axis is positive
proj = x_f @ v
asign = np.sign(proj)
max_ind = np.argmax(np.abs(proj),axis=0)[None,:]
sign = np.take_along_axis(asign,max_ind,axis=0)
proj = proj * sign
return proj
# must be correct https://stackoverflow.com/questions/15230179/how-to-get-the-linear-index-for-a-numpy-array-sub2ind
indices = np.arange(27)
indices3d = indices.reshape((3,3,3))
voxel10 = np.random.normal(size=(3,3,3))
assert voxel10[0,1,2] == voxel10.ravel()[indices3d[0,1,2]]
rot90 = np.rot90
def rotations24(polycube):
# imagine shape is pointing in axis 0 (up)
# 4 rotations about axis 0
yield from rotations4(polycube, 0)
# rotate 180 about axis 1, now shape is pointing down in axis 0
# 4 rotations about axis 0
yield from rotations4(rot90(polycube, 2, axis=1), 0)
# rotate 90 or 270 about axis 1, now shape is pointing in axis 2
# 8 rotations about axis 2
yield from rotations4(rot90(polycube, axis=1), 2)
yield from rotations4(rot90(polycube, -1, axis=1), 2)
# rotate about axis 2, now shape is pointing in axis 1
# 8 rotations about axis 1
yield from rotations4(rot90(polycube, axis=2), 1)
yield from rotations4(rot90(polycube, -1, axis=2), 1)
def rotations4(polycube, axis):
"""List the four rotations of the given cube about the given axis."""
for i in range(4):
yield rot90(polycube, i, axis)
def rot90(m, k=1, axis=2):
"""Rotate an array k*90 degrees in the counter-clockwise direction around the given axis"""
m = np.swapaxes(m, 2, axis)
m = np.rot90(m, k)
m = np.swapaxes(m, 2, axis)
return m
gen = rotations24(indices3d)
res = []
for ind2 in gen:
basis1 = np.array([[i+1,j+1,k+1] for k in range(3) for j in range(3) for i in range(3)]).astype(np.double)
voxel1 = voxel10.copy().reshape(27,1)
voxel2 = voxel1[ind2.reshape(27,)] #np.take(voxel10, ind2).reshape(27,1)
basis1 += 1e-6*np.random.normal(size=(27, 1)) # perturbation
basis2 = basis1
original_diff = np.sum(np.abs(basis1-basis2))
a = pca(voxel1, basis1)
t1 = a[a[:,0].argsort()]
a = pca(voxel2, basis2)
t2 = a[a[:,0].argsort()]
gg= np.abs(t1-t2)
ss= np.sum(np.ravel(gg))
bl= np.all(gg<1e-4)
print('difference before pca %.3f,' % original_diff, 'difference after pca %.3f' % ss,bl)
res.append(bl)
del basis1, basis2
print('correct for %.1f percent of time' % (100*(np.sum(res) / len(res))))
difference before pca 0.000, difference after pca 0.000 True
difference before pca 0.000, difference after pca 0.000 True
difference before pca 0.000, difference after pca 0.000 True
difference before pca 0.000, difference after pca 0.000 True
PS. I want to propose better ordering theme to take into account zero values in the voxel which might confuse previous approach when entire first column of PCA vectors is zero, etc. I propose to sort by vector norms, multiplied by the sign of the sum of elements. Here is tensorflow 2 code:
def infer_shape(x):
x = tf.convert_to_tensor(x)
# If unknown rank, return dynamic shape
if x.shape.dims is None:
return tf.shape(x)
static_shape = x.shape.as_list()
dynamic_shape = tf.shape(x)
ret = []
for i in range(len(static_shape)):
dim = static_shape[i]
if dim is None:
dim = dynamic_shape[i]
ret.append(dim)
return ret
def merge_last_two_dims(tensor):
shape = infer_shape(tensor)
shape[-2] *= shape[-1]
#shape.pop(1)
shape = shape[:-1]
return tf.reshape(tensor, shape)
def pca(inpt_voxel):
patches = tf.extract_volume_patches(inpt_voxel, ksizes=[1,3,3,3,1], strides=[1, 1,1,1, 1], padding="VALID")
features0 = patches[...,tf.newaxis]*basis
# centered basises
basis1_ = tf.ones(shape=tf.shape(patches[...,tf.newaxis]), dtype=tf.float32)*basis
basis1 = basis1_ - tf.math.divide_no_nan(tf.reduce_sum(features0, axis=-2), tf.reduce_sum(patches, axis=-1)[...,None])[:,:,:,:,None,:]
features = patches[...,tf.newaxis]*basis1
features_centered_basis = features - tf.reduce_mean(features, axis=-2)[:,:,:,:,None,:]
x = features_centered_basis
m = tf.cast(x.get_shape()[-2], tf.float32)
cov = tf.matmul(x,x,transpose_a=True)/(m - 1)
e,v = tf.linalg.eigh(cov,name="eigh")
proj = tf.matmul(x,v,transpose_b=False)
asign = tf.sign(proj)
max_ind = tf.argmax(tf.abs(proj),axis=-2)[:,:,:,:,None,:]
sign = tf.gather(asign,indices=max_ind, batch_dims=4, axis=-2)
sign = tf.linalg.diag_part(sign)
proj = proj * sign
# But we can have 1st coordinate zero. In this case,
# other coordinates become ambiguous
#s = tf.argsort(proj[...,0], axis=-1)
# sort by l2 vector norms, multiplied by signs of sums
sum_signs = tf.sign(tf.reduce_sum(proj, axis=-1))
norms = tf.norm(proj, axis=-1)
s = tf.argsort(sum_signs*norms, axis=-1)
proj = tf.gather(proj, s, batch_dims=4, axis=-2)
return merge_last_two_dims(proj)
924
Firstly, your pca function is not correct, it should be
def pca(x):
x -= np.mean(x,axis=0)
cov = np.cov(x.T)
e_values, e_vectors = np.linalg.eig(cov)
order = np.argsort(e_values)[::-1]
v = e_vectors[:,order]
return x @ v
You shouldn't transpose the e_vectors[:,order] because we want each column of the v array is an eigenvector, therefore, x @ v will be projections of x on those eigenvectors.
Secondly, I think you misunderstand the meaning of rotation. It is not voxel1 that should be rotated, but the basis1. If you rotate (by taking transposition) voxel1, what you really do is to rearrange the indices of grid points, while the coordinates of the points basis1 are not changed.
In order to rotate the points (around the z axis for example), you can first define a function to calculate the rotation matrix given an angle
rotate_mat = lambda theta: np.array([[np.cos(theta),-np.sin(theta),0.],[np.sin(theta),np.cos(theta),0.],[0.,0.,1.]])
with the rotation matrix generated by this function, you can rotate the array basis1 to create another array basis2
basis2 = basis1 @ rotate_mat(np.deg2rad(angle))
Now it comes to the title of your question "Why my PCA is not invariant to rotation and axis swap?", from this post, the PCA result is not unique, you can actually run a test to see this
import numpy as np
np.random.seed(10)
rotate_mat = lambda theta: np.array([[np.cos(theta),-np.sin(theta),0.],[np.sin(theta),np.cos(theta),0.],[0.,0.,1.]])
def pca(x):
x -= np.mean(x,axis=0)
cov = np.cov(x.T)
e_values, e_vectors = np.linalg.eig(cov)
order = np.argsort(e_values)[::-1]
v = e_vectors[:,order]
return x @ v
angles = np.linspace(0,90,91)
res = []
for angle in angles:
voxel1 = np.random.normal(size=(3,3,3))
voxel2 = voxel1.copy()
voxel1 = voxel1.reshape(27,1)
voxel2 = voxel2.reshape(27,1)
basis1 = np.array([[i+1,j+1,k+1] for k in range(3) for j in range(3) for i in range(3)])
# basis2 = np.hstack((-basis1[:,1][:,None],basis1[:,0][:,None],-basis1[:,2][:,None]))
basis2 = basis1 @ rotate_mat(np.deg2rad(angle))
voxel1 = voxel1*basis1
voxel2 = voxel2*basis2
print(angle,np.all(np.abs(pca(voxel1) - pca(voxel2) < 1e-6)))
res.append(np.all(np.abs(pca(voxel1) - pca(voxel2) < 1e-6)))
print()
print(np.sum(res) / len(angles))
After you run this script, you will see that in only 21% of times the two PCA results are the same.
UPDATE
I think instead of focusing on the eigenvectors of the principal components, you can instead focus on the projections. For two clouds of points, even though they are essentially the same, the eigenvectors can be drastically different. Therefore, hardcoding in order to somehow let the two sets of eigenvectors to be the same is a very difficult task.
However, based on this post, for the same cloud of points, two sets of eigenvectors can be different only up to a minus sign. Therefore, the projections upon the two sets of eigenvectors are also different only up to a minus sign. This actually offers us an elegant solution, for the projections along an eigenvector (principal axis), all we need to do is to switch the sign of the projections so that the projection with the largest absolute value along that principal axis is positive.
import numpy as np
from numpy.random import randn
#np.random.seed(20)
rotmat_z = lambda theta: np.array([[np.cos(theta),-np.sin(theta),0.],[np.sin(theta),np.cos(theta),0.],[0.,0.,1.]])
rotmat_y = lambda theta: np.array([[np.cos(theta),0.,np.sin(theta)],[0.,1.,0.],[-np.sin(theta),0.,np.cos(theta)]])
rotmat_x = lambda theta: np.array([[1.,0.,0.],[0.,np.cos(theta),-np.sin(theta)],[0.,np.sin(theta),np.cos(theta)]])
# based on https://en.wikipedia.org/wiki/Rotation_matrix
rot_mat = lambda alpha,beta,gamma: rotmat_z(alpha) @ rotmat_y(beta) @ rotmat_x(gamma)
deg2rad = lambda alpha,beta,gamma: [np.deg2rad(alpha),np.deg2rad(beta),np.deg2rad(gamma)]
def pca(feat, x):
# pca with attemt to create convention on sign changes
x_c = x - np.mean(x,axis=0)
x_f = feat * x
x_f -= np.mean(x_f,axis=0)
cov = np.cov(x_f.T)
e_values, e_vectors = np.linalg.eig(cov)
order = np.argsort(e_values)[::-1]
v = e_vectors[:,order]
# here is the solution, we switch the sign of the projections
# so that the projection with the largest absolute value along a principal axis is positive
proj = x_f @ v
asign = np.sign(proj)
max_ind = np.argmax(np.abs(proj),axis=0)[None,:]
sign = np.take_along_axis(asign,max_ind,axis=0)
proj = proj * sign
return proj
ref_angles = np.linspace(0.0,90.0,10)
angles = [[alpha,beta,gamma] for alpha in ref_angles for beta in ref_angles for gamma in ref_angles]
voxel1 = np.random.normal(size=(3,3,3))
res = []
for angle in angles:
voxel2 = voxel1.copy()
voxel1 = voxel1.reshape(27,1)
voxel2 = voxel2.reshape(27,1)
basis1 = np.array([[i+1,j+1,k+1] for k in range(3) for j in range(3) for i in range(3)]).astype(np.double)
basis1 = basis1 + 1e-4 * randn(27,3) # perturbation
basis2 = basis1 @ rot_mat(*deg2rad(*angle))
original_diff = np.sum(np.abs(basis1-basis2))
gg= np.abs(pca(voxel1, basis1) - pca(voxel1, basis2))
ss= np.sum(np.ravel(gg))
bl= np.all(gg<1e-4)
print('difference before pca %.3f,' % original_diff, 'difference after pca %.3f' % ss,bl)
res.append(bl)
del basis1, basis2
print('correct for %.1f percent of time' % (100*(np.sum(res) / len(res))))
As you can see by running this script, the projections on the principal axis are the same, this means we have resolved the issue of PCA results being not unique.
Reply to EDIT 3
As for the new issue you raised, I think you missed an important point, it is the projections of the cloud of points onto the principal axes that are invariant, not anything else. Therefore, if you rotate voxel1 and obtain voxel2, they are the same in the sense that their own respective projections onto the principal axes of the cloud of points are the same, it actually does not make too much sense to compare pca(voxel1,basis1) with pca(voxel2,basis1).
Furthermore, the method rotate of scipy.ndimage actually changes information, as you can see by running this script
image1 = np.linspace(1,100,100).reshape(10,10)
image2 = scipy.ndimage.rotate(image1, 45, mode='nearest', axes=(0, 1), reshape=False)
image3 = scipy.ndimage.rotate(image2, -45, mode='nearest', axes=(0, 1), reshape=False)
fig,ax = plt.subplots(nrows=1,ncols=3,figsize=(12,4))
ax[0].imshow(image1)
ax[1].imshow(image2)
ax[2].imshow(image3)
The output image is
As you can see the matrix after rotation is not the same as the original one, some information of the original matrix is changed.
Reply to EDIT 4
Actually, we are almost there, the two pca results are different because we are comparing pca components for different points.
indices = np.arange(27)
indices3d = indices.reshape((3,3,3))
# apply rotations to the indices, it is not weights yet
gen = rotations24(indices3d)
# construct the weights
voxel10 = np.random.normal(size=(3,3,3))
res = []
count = 0
for ind2 in gen:
count += 1
# ind2 is the array of indices after rotation
# reindex the weights with the indices after rotation
voxel1 = voxel10.copy().reshape(27,1)
voxel2 = voxel1[ind2.reshape(27,)]
# basis1 is the array of coordinates where the points are
basis1 = np.array([[i+1,j+1,k+1] for k in range(3) for j in range(3) for i in range(3)]).astype(np.double)
basis1 += 1e-4*np.random.normal(size=(27, 1))
# reindex the coordinates with the indices after rotation
basis2 = basis1[ind2.reshape(27,)]
# add a slight modification to pca, return the axes also
pca1,v1 = pca(voxel1,basis1)
pca2,v2 = pca(voxel2,basis2)
# sort the principal components before comparing them
pca1 = np.sort(pca1,axis=0)
pca2 = np.sort(pca2,axis=0)
gg= np.abs(pca1 - pca2)
ss= np.sum(np.ravel(gg))
bl= np.all(gg<1e-4)
print('difference after pca %.3f' % ss,bl)
res.append(bl)
del basis1, basis2
print('correct for %.1f percent of time' % (100*(np.sum(res) / len(res))))
Running this script, you will find, for each rotation, the two sets of principal axes are different only up to a minus sign. The two sets of pca results are different because the indices of the cloud of points before and after rotation are different (since you apply rotation to the indices). If you sort the pca results before comparing them, you will find the two pca results are exactly the same.
Summary
The answer to this question can be divided into two parts. In the first part, the rotation is applied to the basis (the coordinates of points), while the indices and the corresponding weights are unchanged. In the second part, the rotation is applied to the indices, then the weights and the basis are rearranged with the new indices. For both of the two parts, the solution pca function is the same
def pca(feat, x):
# pca with attemt to create convention on sign changes
x_c = x - np.mean(x,axis=0)
x_f = feat * x
x_f -= np.mean(x_f,axis=0)
cov = np.cov(x_f.T)
e_values, e_vectors = np.linalg.eig(cov)
order = np.argsort(e_values)[::-1]
v = e_vectors[:,order]
# here is the solution, we switch the sign of the projections
# so that the projection with the largest absolute value along a principal axis is positive
proj = x_f @ v
asign = np.sign(proj)
max_ind = np.argmax(np.abs(proj),axis=0)[None,:]
sign = np.take_along_axis(asign,max_ind,axis=0)
proj = proj * sign
return proj
The idea of this function is, instead of matching the principal axes, we can match the principal components since it is the principal components that are rotationally invariant after all.
Based on this function pca, the first part of this answer is easy to understand, since the indices of the points are unchanged while we only rotate the basis. In order to understand the second part of this answer (Reply to EDIT 5), we must first understand the function rotations24. This function rotates the indices rather than the coordinates of the points, therefore, if we stay at the same position observing the points, we will feel that the positions of the points are changed.
With this in mind, it is not hard to understand Reply to EDIT 5.
Actually, the function pca in this answer can be applied to more general cases, for example (we rotate the indices)
num_of_points_per_dim = 10
num_of_points = num_of_points_per_dim ** 3
indices = np.arange(num_of_points)
indices3d = indices.reshape((num_of_points_per_dim,num_of_points_per_dim,num_of_points_per_dim))
voxel10 = 100*np.random.normal(size=(num_of_points_per_dim,num_of_points_per_dim,num_of_points_per_dim))
gen = rotations24(indices3d)
res = []
for ind2 in gen:
voxel1 = voxel10.copy().reshape(num_of_points,1)
voxel2 = voxel1[ind2.reshape(num_of_points,)]
basis1 = 100*np.random.rand(num_of_points,3)
basis2 = basis1[ind2.reshape(num_of_points,)]
pc1 = np.sort(pca(voxel1, basis1),axis=0)
pc2 = np.sort(pca(voxel2, basis2),axis=0)
gg= np.abs(pc1-pc2)
ss= np.sum(np.ravel(gg))
bl= np.all(gg<1e-4)
print('difference after pca %.3f' % ss,bl)
res.append(bl)
del basis1, basis2
print('correct for %.1f percent of time' % (100*(np.sum(res) / len(res))))
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With
