icosahedron

icosahedron

#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
Icosahedron generation

Imported from vispy (http://vispy.org)
"""

import numpy as np
from PyQt5.QtCore import *
from PyQt5.QtGui import *
from PyQt5.Qt import *
import pyqtgraph as pg
from pyqtgraph.opengl import *

def icosahedron(radius, subdivisions):
    # golden ratio
    t = (1.0 + np.sqrt(5.0))/2.0

    # vertices of a icosahedron
    verts = [(-1, t, 0),
             (1, t, 0),
             (-1, -t, 0),
             (1, -t, 0),
             (0, -1, t),
             (0, 1, t),
             (0, -1, -t),
             (0, 1, -t),
             (t, 0, -1),
             (t, 0, 1),
             (-t, 0, -1),
             (-t, 0, 1)]

    # faces of the icosahedron
    faces = [(0, 11, 5),
             (0, 5, 1),
             (0, 1, 7),
             (0, 7, 10),
             (0, 10, 11),
             (1, 5, 9),
             (5, 11, 4),
             (11, 10, 2),
             (10, 7, 6),
             (7, 1, 8),
             (3, 9, 4),
             (3, 4, 2),
             (3, 2, 6),
             (3, 6, 8),
             (3, 8, 9),
             (4, 9, 5),
             (2, 4, 11),
             (6, 2, 10),
             (8, 6, 7),
             (9, 8, 1)]

    def midpoint(v1, v2):
        return ((v1[0]+v2[0])/2, (v1[1]+v2[1])/2, (v1[2]+v2[2])/2)

    # subdivision
    for _ in range(subdivisions):
        for idx in range(len(faces)):
            i, j, k = faces[idx]
            a, b, c = verts[i], verts[j], verts[k]
            ab, bc, ca = midpoint(a, b), midpoint(b, c), midpoint(c, a)
            verts += [ab, bc, ca]
            ij, jk, ki = len(verts)-3, len(verts)-2, len(verts)-1
            faces.append([i, ij, ki])
            faces.append([ij, j, jk])
            faces.append([ki, jk, k])
            faces[idx] = [jk, ki, ij]
    verts = np.array(verts)
    faces = np.array(faces)
    # make each vertex to lie on the sphere
    lengths = np.sqrt((verts*verts).sum(axis=1))
    verts /= lengths[:, np.newaxis]/radius

    return verts, faces

## Test routine
if __name__ == '__main__':
    verts, faces = icosahedron(100, 1)

    app = QApplication([])
    w = GLViewWidget()
    w.show()
    w.setWindowTitle('pyqtgraph example: GLMeshItem')
    w.setCameraPosition(distance=440)

    colors = np.random.random(size=(faces.shape[0], 4))/2
    mesh = MeshData(vertexes=verts, faces=faces, faceColors=colors)
    n = mesh.faceNormals()
    n /= np.linalg.norm(n[0])

    iso = GLMeshItem(meshdata=mesh, smooth=False)
    #w.addItem(iso)

    votes = np.zeros(faces.shape[0])

    for j in range(300):
        m = np.random.random(3)-np.array([0.5, 0.5, 0.5])
        m /= np.linalg.norm(m)

        line = GLLinePlotItem(pos=np.array([[0,0,0], m*200]), color=(1, 1, 0, 1), width=1)
        #w.addItem(line)

        c = np.zeros(faces.shape[0])
        for i, nn in enumerate(n):
            c[i] = np.dot(nn, m)/np.linalg.norm(nn)

        i = np.argmax(c)
        votes[i] += 1

        line = GLLinePlotItem(pos=np.array([[0,0,0], n[i]*200]), color=(0, 0, 1, 1), width=1)
        #w.addItem(line)

        f = np.array(faces[i])
        v = np.array(verts[f])
        t = GLMeshItem(vertexes=v, faces=np.array([[0, 1, 2]]), faceColors=np.array([np.random.random(4)]), smooth=False)
        w.addItem(t)

    print(votes)

    QApplication.instance().exec_()