python - 使用 PyOpenGL 可视化 3D NumPy 数组

Question

我想创建一个 PyOpenGL/QtOpenGL 小部件，它允许我可视化任意 NumPy 3D 矩阵，与以下 Hinton 图不同，设想为“立方体的立方体”而不是“正方形的正方形”:

虽然我在使用 OpenGL 时遇到了一些困难。到目前为止，这是我的代码:

from OpenGL.GL import *
from OpenGL.GLUT import *
from PyQt4 import QtGui, QtOpenGL
import numpy as np

action_keymap = {
  # 'a': lambda: glTranslate(-1, 0, 0),
  # 'd': lambda: glTranslate( 1, 0, 0),
  # 'w': lambda: glTranslate( 0, 1, 0),
  # 's': lambda: glTranslate( 0,-1, 0),

  'a': lambda: glRotate(-5, 0, 1, 0),
  'd': lambda: glRotate( 5, 0, 1, 0),
  # 'W': lambda: glRotate(-5, 1, 0, 0),
  # 'S': lambda: glRotate( 5, 1, 0, 0),
}

ARRAY = np.ones([3,3,3])

class GLWidget(QtOpenGL.QGLWidget):

  def paintGL(self):
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

    for idx, value in np.ndenumerate(ARRAY):
      rel_pos = np.array(idx)/np.max(ARRAY.shape)
      glTranslate(* rel_pos)
      glutSolidCube(0.9/np.max(ARRAY.shape))
      glTranslate(*-rel_pos)

  def resizeGL(self, w, h):
    glLoadIdentity()
    glRotate(35,1,0,0)
    glRotate(45,0,1,0)

  def initializeGL(self):
    glClearColor(0.1, 0.1, 0.3, 1.0)

  def keyPressEvent(self, event):
    action = action_keymap.get(str(event.text()))
    if action:
        action()
    self.updateGL()

  def mousePressEvent(self, event):
    super().mousePressEvent(event)
    self.press_point = event.pos()

  def mouseMoveEvent(self, event):
    super().mouseMoveEvent(event)
    motion = event.pos()-self.press_point
    self.press_point = event.pos()
    glRotate(motion.x(),0,1,0)
    glRotate(motion.y(),1,0,0)
    self.updateGL()

if __name__ == '__main__':
  app = QtGui.QApplication(sys.argv)

  w = GLWidget()
  w.show()

  sys.exit(app.exec_())

我的问题如下:

1) 照明。我一直在阅读有关照明和 Material 的信息，但我似乎无法在某处获得简单的灯光来使形状更加清晰。我想要最简单、最基本的可能光来区分正方形，而不是将它们的所有面都呈现为纯白色。我知道如何更改颜色，但这并不能缓解问题。 我可以用什么最简单的光照射这个格子来让子组件更清晰？

2) 它很慢。我将计算数学以实现正方形的正确定位和调整大小，但我想知道是否有一种方法可以矢量化该过程(毕竟，它只是将索引转换为平移并将值转换为立方体数组中每个元素的大小)。我应该在 cpp 中编写扩展，用 ctypes 包装我的代码，还是有办法明确地将工作外包给 OpenGL？ 从 Python 向 OpenGL 发送重复性任务的标准方法是什么？

Answer 1

此任务非常适合 Instancing .通过实例化，一个对象可以被渲染多次。

在这种情况下，实例化用于为 3d 的每个元素渲染立方体 NumPy数组。

假设我们有以下随机值在 [0, 1] 范围内的 3D 数组 (array3d):

shape = [5, 4, 6]
number_of = shape[0] * shape[1] * shape[2]  
array3d = np.array(np.random.rand(number_of), dtype=np.float32).reshape(shape)

对于数组的每个元素，必须渲染一个网格(立方体)实例:

例如

number_of = array3d.shape[0] * array3d.shape[1] * array3d.shape[2]  
glDrawElementsInstanced(GL_TRIANGLES, self.__no_indices, GL_UNSIGNED_INT, None, number_of)

数组可以加载到 3D 纹理(glTexImage3D):

glActiveTexture(GL_TEXTURE1)
tex3DObj = glGenTextures(1)
glBindTexture(GL_TEXTURE_3D, tex3DObj)
glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAX_LEVEL, 0)
glTexImage3D(GL_TEXTURE_3D, 0, GL_R16F, *array3d.shape, 0, GL_RED, GL_FLOAT, array3d)

在单个立方体的顶点着色器中，可以通过 3D 纹理的维度(等于 3D 数组的形状)和 gl_InstanceID 计算实例变换矩阵。元素立方体。
元素立方体根据 3D 纹理中元素的值进一步缩放。

假设顶点着色器具有 §D 纹理采样器统一 u_array3D 和顶点坐标属性 a_pos:

in vec3 a_pos;
uniform sampler3D u_array3D;

纹理的维度可以通过textureSize得到:

ivec3 dim = textureSize(u_array3D, 0);

与维度和 gl_InstanceID ，可以计算元素的索引:

ivec3 inx = ivec3(0);
inx.z = gl_InstanceID / (dim.x * dim.y);
inx.y = (gl_InstanceID - inx.z * dim.x * dim.y) / dim.x;
inx.x = gl_InstanceID - inx.z * dim.x * dim.y - inx.y * dim.x;

并且可以获取元素的值(texelFetch):

float value = texelFetch(u_array3D, inx, 0).x;

最后可以计算出依赖于元素索引和元素值的实例变换矩阵:

vec3 scale = 1.0 / vec3(dim);
scale = vec3(min(scale.x, min(scale.y, scale.z)));
vec3 trans = 2 * scale * (vec3(inx) - vec3(dim-1) / 2.0);
mat4 instanceMat = mat4(
    vec4(scale.x * cube_scale, 0.0, 0.0, 0.0),
    vec4(0.0, scale.y * cube_scale, 0.0, 0.0),
    vec4(0.0, 0.0, scale.z * cube_scale, 0.0),
    vec4(trans, 1.0)
);

vec4 instance_pos = instanceMat * vec4(a_pos, 1.0);

该值还可以通过立方体的颜色进行可视化。为此，范围 [0.0, 1.0] 中的浮点值被转换为 HSV 中的 RGB 颜色。颜色范围:

vec3 HUEtoRGB(in float H)
{
    float R = abs(H * 6.0 - 3.0) - 1.0;
    float G = 2.0 - abs(H * 6.0 - 2.0);
    float B = 2.0 - abs(H * 6.0 - 4.0);
    return clamp( vec3(R,G,B), 0.0, 1.0 );
}

vec3 color = HUEtoRGB(0.66 * (1-0 - value));

---

另见 OpenGL - Python examples

纯NumPy/PyOpenGL示例程序。数组的值随机变化:

import numpy as np
from OpenGL.GLUT import *
from OpenGL.GL import *
from OpenGL.GL.shaders import *

class MyWindow:

    __glsl_vert = """
        #version 450 core

        layout (location = 0) in vec3 a_pos;
        layout (location = 1) in vec3 a_nv;
        layout (location = 2) in vec4 a_col;

        out vec3 v_pos;
        out vec3 v_nv;
        out vec4 v_color;

        layout (binding = 1) uniform sampler3D u_array3D;

        uniform mat4 u_proj; 
        uniform mat4 u_view; 
        uniform mat4 u_model; 

        vec3 HUEtoRGB(in float H)
        {
            float R = abs(H * 6.0 - 3.0) - 1.0;
            float G = 2.0 - abs(H * 6.0 - 2.0);
            float B = 2.0 - abs(H * 6.0 - 4.0);
            return clamp( vec3(R,G,B), 0.0, 1.0 );
        }

        void main()
        {
            ivec3 dim = textureSize(u_array3D, 0);

            vec3 scale = 1.0 / vec3(dim);
            scale = vec3(min(scale.x, min(scale.y, scale.z)));
            
            ivec3 inx = ivec3(0);
            inx.z = gl_InstanceID / (dim.x * dim.y);
            inx.y = (gl_InstanceID - inx.z * dim.x * dim.y) / dim.x;
            inx.x = gl_InstanceID - inx.z * dim.x * dim.y - inx.y * dim.x;
            float value = texelFetch(u_array3D, inx, 0).x;

            vec3 trans = 2 * scale * (vec3(inx) - vec3(dim-1) / 2.0);
            mat4 instanceMat = mat4(
                vec4(scale.x * value, 0.0, 0.0, 0.0),
                vec4(0.0, scale.y * value, 0.0, 0.0),
                vec4(0.0, 0.0, scale.z * value, 0.0),
                vec4(trans, 1.0)
            );

            mat4 model_view = u_view * u_model * instanceMat;
            mat3 normal     = transpose(inverse(mat3(model_view)));
            
            vec4 view_pos   = model_view * vec4(a_pos.xyz, 1.0);

            v_pos       = view_pos.xyz;
            v_nv        = normal * a_nv;  
            v_color     = vec4(HUEtoRGB(0.66 * (1-0 - value)), 1.0);
            gl_Position = u_proj * view_pos;
        }
    """

    __glsl_frag = """
        #version 450 core
        
        out vec4 frag_color;
        in  vec3 v_pos;
        in  vec3 v_nv;
        in  vec4 v_color;

        void main()
        {
            vec3  N    = normalize(v_nv);
            vec3  V    = -normalize(v_pos);
            float ka   = 0.1;
            float kd   = max(0.0, dot(N, V)) * 0.9;
            frag_color = vec4(v_color.rgb * (ka + kd), v_color.a);
        }
    """

    def __init__(self, w, h):
        
        self.__caption = 'OpenGL Window'
        self.__vp_valid = False
        self.__vp_size = [w, h]

        glutInit()
        glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH)
        glutInitWindowSize(self.__vp_size[0], self.__vp_size[1])
        self.__glut_wnd = glutCreateWindow(self.__caption)

        self.__program = compileProgram( 
            compileShader( self.__glsl_vert, GL_VERTEX_SHADER ),
            compileShader( self.__glsl_frag, GL_FRAGMENT_SHADER ),
        )
        self.___attrib = { a : glGetAttribLocation (self.__program, a) for a in ['a_pos', 'a_nv', 'a_col'] }
        print(self.___attrib)
        self.___uniform = { u : glGetUniformLocation (self.__program, u) for u in ['u_model', 'u_view', 'u_proj'] }
        print(self.___uniform)

        v = [[-1,-1,1], [1,-1,1], [1,1,1], [-1,1,1], [-1,-1,-1], [1,-1,-1], [1,1,-1], [-1,1,-1]]
        c = [[1.0, 0.0, 0.0], [1.0, 0.5, 0.0], [1.0, 0.0, 1.0], [1.0, 1.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]
        n = [[0,0,1], [1,0,0], [0,0,-1], [-1,0,0], [0,1,0], [0,-1,0]]
        e = [[0,1,2,3], [1,5,6,2], [5,4,7,6], [4,0,3,7], [3,2,6,7], [1,0,4,5]]
        index_array = [si*4+[0, 1, 2, 0, 2, 3][vi] for si in range(6) for vi in range(6)]
        attr_array = []
        for si in range(len(e)):
            for vi in e[si]:
                attr_array += [*v[vi], *n[si], *c[si], 1]
        
        self.__no_vert = len(attr_array) // 10
        self.__no_indices = len(index_array)
        vertex_attributes = np.array(attr_array, dtype=np.float32)
        indices = np.array(index_array, dtype=np.uint32)
        
        self.__vao = glGenVertexArrays(1)
        self.__vbo, self.__ibo = glGenBuffers(2)
        
        glBindVertexArray(self.__vao)

        glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, self.__ibo)
        glBufferData(GL_ELEMENT_ARRAY_BUFFER, indices, GL_STATIC_DRAW)

        glBindBuffer(GL_ARRAY_BUFFER, self.__vbo)
        glBufferData(GL_ARRAY_BUFFER, vertex_attributes, GL_STATIC_DRAW)

        float_size = vertex_attributes.itemsize  
        glVertexAttribPointer(0, 3, GL_FLOAT, False, 10*float_size, None)
        glVertexAttribPointer(1, 3, GL_FLOAT, False, 10*float_size, c_void_p(3*float_size))
        glVertexAttribPointer(2, 4, GL_FLOAT, False, 10*float_size, c_void_p(6*float_size))
        glEnableVertexAttribArray(0)
        glEnableVertexAttribArray(1)
        glEnableVertexAttribArray(2)

        glEnable(GL_DEPTH_TEST)
        glUseProgram(self.__program)

        shape = [5, 4, 6]
        number_of = shape[0] * shape[1] * shape[2]  
        self.array3d = np.array(np.random.rand(number_of), dtype=np.float32).reshape(shape)

        glActiveTexture(GL_TEXTURE1)
        self.tex3DObj = glGenTextures(1)
        glBindTexture(GL_TEXTURE_3D, self.tex3DObj)
        glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAX_LEVEL, 0)
        glTexImage3D(GL_TEXTURE_3D, 0, GL_R16F, *self.array3d.shape, 0, GL_RED, GL_FLOAT, self.array3d)

        glutReshapeFunc(self.__reshape)
        glutDisplayFunc(self.__mainloop)

    def run(self):
        self.__starttime = 0
        self.__starttime = self.elapsed_ms()
        glutMainLoop()

    def elapsed_ms(self):
    return glutGet(GLUT_ELAPSED_TIME) - self.__starttime

    def __reshape(self, w, h):
        self.__vp_valid = False

    def __mainloop(self):

        number_of = self.array3d.shape[0] * self.array3d.shape[1] * self.array3d.shape[2]  
        rand = (np.random.rand(number_of) - 0.5) * 0.05
        self.array3d = np.clip(np.add(self.array3d, rand.reshape(self.array3d.shape)), 0, 1)
        glTexSubImage3D(GL_TEXTURE_3D, 0, 0, 0, 0, *self.array3d.shape, GL_RED, GL_FLOAT, self.array3d)

        if not self.__vp_valid:
            self.__vp_size = [glutGet(GLUT_WINDOW_WIDTH), glutGet(GLUT_WINDOW_HEIGHT)]
            self.__vp_valid = True
            glViewport(0, 0, self.__vp_size[0], self.__vp_size[1])

        aspect, ta, near, far = self.__vp_size[0]/self.__vp_size[1], np.tan(np.radians(90.0) / 2), 0.1, 10
        proj = np.array(((1/ta/aspect, 0, 0, 0), (0, 1/ta, 0, 0), (0, 0, -(far+near)/(far-near), -1), (0, 0, -2*far*near/(far-near), 0)), np.float32)
        
        view = np.array(((1, 0, 0, 0), (0, 0, -1, 0), (0, 1, 0, 0), (0, 0, -2, 1)), np.float32)
        c, s = (f(np.radians(30.0)) for f in [np.cos, np.sin])
        viewRotX = np.array(((1, 0, 0, 0), (0, c, s, 0), (0, -s, c, 0), (0, 0, 0, 1)), np.float32)
        view = np.matmul(viewRotX, view)
        
        c1, s1, c2, s2, c3, s3 = (f(self.elapsed_ms() * np.pi * 2 / tf) for tf in [5000.0, 7333.0, 10000.0] for f in [np.cos, np.sin])
        rotMatZ = np.array(((c3, s3, 0, 0), (-s3, c3, 0, 0), (0, 0, 1, 0), (0, 0, 0, 1)), np.float32)
        model = rotMatZ
    
        glUniformMatrix4fv(self.___uniform['u_proj'], 1, GL_FALSE, proj )
        glUniformMatrix4fv(self.___uniform['u_view'], 1, GL_FALSE, view )
        glUniformMatrix4fv(self.___uniform['u_model'], 1, GL_FALSE, model )

        glClearColor(0.2, 0.3, 0.3, 1.0)
        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
        
        glDrawElementsInstanced(GL_TRIANGLES, self.__no_indices, GL_UNSIGNED_INT, None, number_of)

        glutSwapBuffers()
        glutPostRedisplay()

window = MyWindow(800, 600)
window.run()