How can an almost arbitrary plane in a 3D dataset be plotted by matplotlib?

Question

There is an array containing 3D data of shape e.g. (64,64,64), how do you plot a plane given by a point and a normal (similar to hkl planes in crystallography), through this dataset? Similar to what can be done in MayaVi by rotating a plane through the data.

The resulting plot will contain non-square planes in most cases. Can those be done with matplotlib (some sort of non-rectangular patch)?

Edit: I almost solved this myself (see below) but still wonder how non-rectangular patches can be plotted in matplotlib...?

Edit: Due to discussions below I restated the question.

Answer 1

This is funny, a similar question I replied to just today. The way to go is: interpolation. You can use griddata from scipy.interpolate:

Griddata

This page features a very nice example, and the signature of the function is really close to your data.

You still have to somehow define the points on you plane for which you want to interpolate the data. I will have a look at this, my linear algebra lessons where a couple of years ago

Answer 2

I have the penultimate solution for this problem. Partially solved by using the second answer to Plot a plane based on a normal vector and a point in Matlab or matplotlib :

# coding: utf-8
import numpy as np
from matplotlib.pyplot import imshow,show

A=np.empty((64,64,64)) #This is the data array
def f(x,y):
    return np.sin(x/(2*np.pi))+np.cos(y/(2*np.pi))
xx,yy= np.meshgrid(range(64), range(64))
for x in range(64):
    A[:,:,x]=f(xx,yy)*np.cos(x/np.pi)

N=np.zeros((64,64)) 
"""This is the plane we cut from A. 
It should be larger than 64, due to diagonal planes being larger. 
Will be fixed."""

normal=np.array([-1,-1,1]) #Define cut plane here. Normal vector components restricted to integers
point=np.array([0,0,0])
d = -np.sum(point*normal)

def plane(x,y): # Get plane's z values
    return (-normal[0]*x-normal[1]*y-d)/normal[2]

def getZZ(x,y): #Get z for all values x,y. If z>64 it's out of range
    for i in x:
        for j in y:
            if plane(i,j)<64:
                N[i,j]=A[i,j,plane(i,j)]

getZZ(range(64),range(64))
imshow(N, interpolation="Nearest")
show()

It's not the ultimate solution since the plot is not restricted to points having a z value, planes larger than 64 * 64 are not accounted for and the planes have to be defined at (0,0,0).