Well basically title is what I'm trying to do except I would also like to
sum up (integrate) the projected values in bins.
In other words it could be described as a 2D heatmap of a 3D point set.
In another other words, imagine a normal distribution projected from the z
axis to the xy plane. What you see is basically a 2D normal distribution.
Moving an azimuthal angle theta (like in spherical coordinates) away from
the z axis when you project the same normal distribution what you would get
is not the "normal" normal dist but a skewed normal distribution.
I thought I can get the same results by rotating the normal distribution,
projecting it on xy plane and then doing a heatmap, hexbin, histogram or
something similar. Here's the rotation an projection via the contourf
method, but I can't reproduce the projections myself.
I create an x, y and z as 2D arrays (xy meshgird of coordinates, and z as a
Cartesian product of 1D normal dist.). I rotate them by multiplying the
rotation matrix with the column vector of every point on the grid just like
in basic lin alg. No issues.
Rotation matrix is created in the rotation_matrix_weave but for the people
that have too modern scipy's/mpl's the rotation_matrix_numpy does the same
thing it just takes more time.
The points are rotated in the RotateFunc (where the meshgrid and Cartesian
product is done too).
I plot the rotated normal dist. and use the contourf functions to show the
projections on xy and yz planes. The kind of projection that is visible in
the xy plane is the one I'd like to get. In all fairness because these are
projections of points if you would just rotate an empty xy plane what you
would get as a projection back on to the xy plane is a gradient. So I
rotate an "empty" plane and subtract it from the rotated normal projection
to get rid of the gradient and just keep the skewed gaussian.
But this is where all the problems start. I can't figure out how to get
that projection. I tried finding where the countourf functions are defined,
so I could try to figure out how it does it but couldn't make out the code
(it's fairly big) and it doesn't really do what I need it to.
Easier approach is to just construct the projections over pcolor or
pcolormesh (much faster) as I've added at the bottom of the script. But the
online wisdom on how to reconstruct the arrays won't work because
pcolormesh/pcolor understands that "indices aren't coordinates". That is to
say, pcolormesh doesn't simply just plot rz, but it plots the value rz[i,j]
at the coordinates ( rx[i,j], ry[i,j] ) in this new 2D array that is the
final image.
I've been fishing for help online and in IRC for 2-3 days now and I've been
lost for over a week, by this point honestly, whatever works. I've resorted
to saving the graph as an image and reading that in with OpenCv. But the
problem is they are all "aspected" so from the original (N, N) array you
end up with a (1000, 800) image and I really don't like the possibility
that the original data can change like that. I mean contourf, pcolor and
pcolormesh can do it, there's got to be a way to extract the 2D array or at
least some kind of bin-value pairs or something from them.
So to reiterate:
1) can someone show me how to get the end graph shown by pcolormesh as an
array?
2) if that's not possible is pcolor any different?
3) if it's not, does anyone have insight on how pcolormesh/pcolor can end
up drawing the graph correctly and can I get pointers on how to do what
they do myself?
4) is there perhaps some other method/function I'm not looking at (i.e
hexbin) that can be used to get the result?
5) if that's not possible at all, could it be possible to change from the
current representation to i.e. [ [x, y, z]....] of the points and use that
to get a histogram/hexbin? Those are the xyz coordinates of points, that is
[ [ rx[0,0] , ry[0,0] , rz[0,0] ], ...... [ rx[i,j] , ry[i,j] , rz[i,j]
], ..... [ rx[-1,-1] , ry[-1,-1] , rz[-1,-1] ] ].
6) am I stupid?
Sorry for the long mail, thanks for all the help.
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy import stats
from matplotlib import cm
from math import cos, sin, sqrt
from scipy import weave
def rotation_matrix_numpy(axis, theta):
mat = np.eye(3,3)
axis = axis/sqrt(np.dot(axis, axis))
a = cos(theta/2.)
b, c, d = -axis*sin(theta/2.)
return np.array([[a*a+b*b-c*c-d*d, 2*(b*c-a*d), 2*(b*d+a*c)],
[2*(b*c+a*d), a*a+c*c-b*b-d*d, 2*(c*d-a*b)],
[2*(b*d-a*c), 2*(c*d+a*b), a*a+d*d-b*b-c*c]])
def rotation_matrix_weave(theta, axiss,
mat = [1., 0., 0., 0., 1., 0., 0., 0., 1.]):
support = "#include <math.h>"
code = """
double axis [3] = {axiss[0], axiss[1], axiss[2]};
double x = sqrt(axis[0] * axis[0] + axis[1] * axis[1] + axis[2] *
axis[2]);
double a = cos(theta / 2.0);
double b = -(axis[0] / x) * sin(theta / 2.0);
double c = -(axis[1] / x) * sin(theta / 2.0);
double d = -(axis[2] / x) * sin(theta / 2.0);
mat[0] = a*a + b*b - c*c - d*d;
mat[1] = 2 * (b*c - a*d);
mat[2] = 2 * (b*d + a*c);
mat[3*1 + 0] = 2*(b*c+a*d);
mat[3*1 + 1] = a*a+c*c-b*b-d*d;
mat[3*1 + 2] = 2*(c*d-a*b);
mat[3*2 + 0] = 2*(b*d-a*c);
mat[3*2 + 1] = 2*(c*d+a*b);
mat[3*2 + 2] = a*a+d*d-b*b-c*c;
"""
weave.inline(code, arg_names=["mat", "axiss", "theta"],
support_code = support)
return np.reshape(mat, (3,3))
def RotateFunc(X, Y, Z, theta, axis=[1.0, 0.0, 0.0]):
x, y = np.meshgrid(X, Y)
z = np.outer(Z, Z)
R = rotation_matrix_weave(theta, axis)
rx, ry, rz = np.zeros(x.shape), np.zeros(y.shape), np.zeros(z.shape)
for i in range(x.shape[0]):
for j in range(x.shape[1]):
tmpx, tmpy, tmpz = np.matmul(R, [x[i, j], y[i, j], z[i, j]])
rx[i,j]+=tmpx
ry[i,j]+=tmpy
rz[i,j]+=tmpz
return rx, ry, rz
···
###############################################
### The plot
###############################################
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
xmin, xmax = -5, 5
ymin, ymax = -5, 5
theta = 0.3
N = 1000
x = np.linspace(xmin, xmax, N)
y = np.linspace(ymin, ymax, N)
z = 5.*stats.multivariate_normal().pdf(x)
xp = np.linspace(xmin, xmax, N)
yp = np.linspace(xmin, xmax, N)
zp = np.zeros(x.shape)
rx, ry, rz = RotateFunc(x, y, z, theta)
rx2, ry2, rz2 = RotateFunc(xp, yp, zp, theta)
ax.plot_surface(rx, ry, rz, color='red',alpha=0.65, linewidth=1)
#cset = ax.contourf(rx[:500], ry[:500], rz[:500], zdir='x', offset=-4,
cmap=cm.coolwarm)
cset = ax.contourf(rx, ry-ry2, rz, zdir='y', offset=-4, cmap=cm.coolwarm)
cset = ax.contourf(rx, ry, rz-rz2, zdir='z', offset=-4, cmap=cm.coolwarm)
plt.show()
plt.pcolormesh(rx, ry, rz-rz2)
plt.show()
plt.pcolormesh(rx, rz, ry-ry2)
plt.show()
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/matplotlib-users/attachments/20170125/e4bd2b9d/attachment-0001.html>