Orthogonal projection of 3D point set onto an xy plane

Previous Topic Next Topic
classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
Report Content as Inappropriate

Orthogonal projection of 3D point set onto an xy plane

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]])
    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.pcolormesh(rx, ry, rz-rz2)
plt.pcolormesh(rx, rz, ry-ry2)


Matplotlib-users mailing list
[hidden email]