動作環境

GeForce GTX 1070 (8GB)
ASRock Z170M Pro4S [Intel Z170chipset]
Ubuntu 16.04.4 LTS desktop amd64
TensorFlow v1.7.0
cuDNN v5.1 for Linux
CUDA v8.0
Python 3.5.2
IPython 6.0.0 -- An enhanced Interactive Python.
gcc (Ubuntu 5.4.0-6ubuntu1~16.04.4) 5.4.0 20160609
GNU bash, version 4.3.48(1)-release (x86_64-pc-linux-gnu)
scipy v0.19.1
geopandas v0.3.0
MATLAB R2017b (Home Edition)
ADDA v.1.3b6
gnustep-gui-runtime v0.24.0-3.1
PyMieScatt v1.7.0

概要

SIRIS:
- SIRIS | Geometry > Gaussian random sphere (@ SIRIS) > Output > ,and which points should be connected to form the trinangulation
- SIRIS > VTK > ASCII VTK Files > POLYGONS 1152 4608の4608 > the number of cell list size
SIRISというパッケージにておいてGaussian Random Sphereという形状を生成できる
Gaussian Random Sphereを可視化したい
- POV-Rayでは可視化(dispole表現した物体)は実装済
  - ADDA > tool > make_pov_png_180121.py > v0.1-v0.7 > to make [Povray file] and [png files] of the particle shape from the [IntField-Y] in terms of (slant, top, front, and side) views
グルグルまわしたい
- GitHub STLレンダリングを使う
  - 正四角錐では実装済
    - 可視化 (3D) > 正四角錐 (格子状のdipoleで構成) > GitHub + STLレンダリングでの表示 > make_cubeGroup_180520.py v0.1

今回

SIRISの出力vtk.outからTrianglesの座標情報を読込む
読込んだTrianglesの座標情報を元に、格子状のdipole配置で表現 (以下、**「dipole表現」**と記載)
GitHubに上げてぐりぐり

処理の流れ

vtk.outを読込み、dipole表現する
- volfil_vtkShape_180526.py
dipole表現に対して、STLファイル形式にする
- make_cubeGroup_180520.py
  - 入力: dipole_180526.res
  - 出力: cubeGroup_180527_t0805.stl
GitHubにpush

実装

VTKReader_180526.py v0.2

VTK形式ファイルを読む

VTKReader_180526.py

import numpy as np
import sys

'''
v0.2 May, 27, 2018
  - fix bug: get_triangle_coordinate() returned [number of vertices] mistakenly
  - rename split_to_triangles() to get_triangle_coordinate()
v0.1 May, 26, 2018
  - branched
=== branched from [test_read_vtkout_180526.py] ===
v0.2 May, 26, 2018
  - add Test_main()
  - add get_gaussian_sphere_vertices()
  - add read_polygons_vertices()
  - add split_to_triangles()
  - read polygonIdxs[]
  - read vertices[]
  - add get_polygons_number()
v0.1 May, 26, 2018
  - add get_vertices_number()
'''


def get_vertices_number(instr):
    # instr: e.g. 'POINTS    578 float'
    return int(instr.split()[1])


def get_polygons_number(instr):
    # instr: e.g. 'POLYGONS    1152   4608'
    return int(instr.split()[1])


def get_triangle_coordinate(polygonIndices, vertices):
    ret = []
    for elem in polygonIndices:
        # get vertex indices from elem
        # e.g. "           3           2           6           7"
        # where the 1st [3] is the number of polygon vertices
        ret += [[int(elem[1]), int(elem[2]), int(elem[3])]]
    return ret


def read_polygons_vertices(infile):
    # key to find the corresponding line
    KEY_VERTICES = 'POINTS'
    KEY_POLYGONS = 'POLYGONS'

    dat = np.genfromtxt(infile, dtype='str', delimiter=',')

    cdwnVert = -1  # countdown to read [vertices], >0: read
    vertices = []  # 2d array having (x,y,z) for each
    cdwnPoly = -1  # countdown to read [Polygons], >0: read
    polygonIdxs = []
    for elem in dat:
        if cdwnVert > 0:
            cdwnVert -= 1
            vertices += [elem.split()]
        if cdwnPoly > 0:
            cdwnPoly -= 1
            polygonIdxs += [elem.split()]
        if KEY_VERTICES in elem:
            numVert = get_vertices_number(elem)
            cdwnVert = numVert
        if KEY_POLYGONS in elem:
            numPlygn = get_polygons_number(elem)
            cdwnPoly = numPlygn
    return polygonIdxs, vertices


def get_gaussian_sphere_vertices(infile):
    plygnIdx, verts = read_polygons_vertices(infile)

    tris = get_triangle_coordinate(plygnIdx, verts)

    ret = []
    for atri in tris:
        vrt1 = [float(x) for x in verts[atri[0]]]
        vrt2 = [float(x) for x in verts[atri[1]]]
        vrt3 = [float(x) for x in verts[atri[2]]]
        ret += [np.array([vrt1, vrt2, vrt3])]
    return np.array(ret)


def Test_main():
    IN_FILE = 'vtk.out'  # output of [G-sphere]

    vrts = get_gaussian_sphere_vertices(IN_FILE)
    print(vrts[0])
    print(vrts[-1])


if __name__ == '__main__':
    Test_main()

volfil_vtkShape_180526.py v0.2

vtk.outをもとにdipole表現する
dpl = get_dipole_positions_in_rectangle(expand=2.0, numdiv=20)
- expand=2.0: 対象領域を[-2.0, 2.0]とする (2.0を含む)
- numdiv=20: 対象領域を20x20x20のdipoleで表現する
- expand, numdivは対象形状、分解能に合わせて変更して使う

volfil_vtkShape_180526.py

import numpy as np
import geometry_intersection_3D_180513 as GI3
import VTKReader_180526 as VTKRD
import sys
import os


'''
v0.2 May, 27, 2018
  - get_dipole_positions_in_rectangle*) takes [numdiv] arg
v0.1 May, 26, 2018
  - import [VTKReader_180526.py]
=== branched from [volfil_tetrahedron_180519.py] ===
v0.2 May, 20, 2018
  - output as file for reading with Numpy genfromtxt()
    + add [OUTPUT]
    + remove [OUTPUT] before file save
v0.1 May, 19, 2018
  - add check_vertex_intersecting_ray()
  - import [geometry_intersection_3D_180513]
  - add get_dipole_positions_in_rectangle()
  - add get_regular_tetrahedron_vertices()
'''

OUTFILE = 'dipole_180526.res'  # for Numpy genfromtxt()


def get_regular_tetrahedron_vertices(expand=1.0):
    vrt1 = np.array([expand, expand, expand])  # vertex
    vrt2 = np.array([-expand, -expand, expand])
    vrt3 = np.array([-expand, expand, -expand])
    vrt4 = np.array([expand, -expand, -expand])

    ret = []
    ret += [np.array([vrt1, vrt2, vrt3])]
    ret += [np.array([vrt1, vrt3, vrt4])]
    ret += [np.array([vrt2, vrt3, vrt4])]
    ret += [np.array([vrt1, vrt2, vrt4])]

    return np.array(ret)


def get_dipole_positions_in_rectangle(expand=5.0, numdiv=9):
    mx = np.linspace(-expand, expand, numdiv)
    my = np.linspace(-expand, expand, numdiv)
    mz = np.linspace(-expand, expand, numdiv)
    ret = []
    for idx in range(len(mx)):
        for idy in range(len(my)):
            for idz in range(len(mz)):
                ret += [np.array([mx[idx], my[idy], mz[idz]])]
    return np.array(ret)


def check_vertex_intersecting_ray(start, end, vertices):
    ray = np.array([start, end])
    for vertex in vertices:
        res, retI = GI3.check_intersection(ray, vertex)
        isIntsct = GI3.check_reachingIntersection(ray, vertex)
        if isIntsct:
            return False
    return True


IN_FILE = 'vtk.out'
vtkshp = VTKRD.get_gaussian_sphere_vertices(IN_FILE)
dpl = get_dipole_positions_in_rectangle(expand=2.0, numdiv=20)

startary = np.array([0, 0, 0])

os.remove(OUTFILE)
with open(OUTFILE, 'ab') as wfp:
    for elem in dpl:
        res = check_vertex_intersecting_ray(startary, elem, vtkshp)
        if res:
            print("[%f, %f, %f]," % (elem[0], elem[1], elem[2]))
            wrk = [elem]  # to avoid element-wise-LF for 1-D array
            np.savetxt(wfp, wrk, fmt='%.7f', delimiter=' ')

make_cubeGroup_180520.py v0.2

dipole表現を元にSTLファイルを生成する

make_cubeGroup_180520.py

import numpy as np
import sys

import STLWriter as STLWR
# STLWriter
# code from
# http://code.activestate.com/recipes/578246/
# by Manfred Moitzi
# (renamed as STLWriter.py from recipe-578246-1.py)

#

'''
v0.2 May, 26, 2018
  - change [OUT_FILE]
v0.1 May, 20, 2018
  - read dipole positions from [IN_FILE]
  - add python_list_add()
  - get_cube() takes [size] arg
  - add make_cubeGroup()
'''


def python_list_add(in1, in2):
    wrk = np.array(in1) + np.array(in2)
    return wrk.tolist()


def make_cubeGroup():
    def get_cube(size=3, origin=[0, 0, 0]):
        # cube corner points
        s = size
        p1 = python_list_add(origin, (0, 0, 0))
        p2 = python_list_add(origin, (0, 0, s))
        p3 = python_list_add(origin, (0, s, 0))
        p4 = python_list_add(origin, (0, s, s))
        p5 = python_list_add(origin, (s, 0, 0))
        p6 = python_list_add(origin, (s, 0, s))
        p7 = python_list_add(origin, (s, s, 0))
        p8 = python_list_add(origin, (s, s, s))

        print(p1)

        # define the 6 cube faces
        # faces just lists of 3 or 4 vertices
        #
        # [NOTE]
        # vertices must be defined [anticlockwise] to be seen on GitHub
        # otherwise, the face is not drawn
        #
        return [
            [p1, p3, p7, p5],  # bottom
            [p1, p5, p6, p2],
            [p5, p7, p8, p6],
            [p7, p3, p4, p8],  # rear
            [p1, p2, p4, p3],  # left
            [p2, p6, p8, p4],
        ]

    # 1. Regular Tetrahedra
    # IN_FILE = 'dipole_180520.res'  # from [volfil_tetrahedron_180519.py] v0.2
    # OUT_FILE = 'cubeGroup_180520_t0950.stl'
    # DIPOLE_SIZE = 5

    # 2. VTK from Gaussian Sphere code
    IN_FILE = 'dipole_180526.res'  # from [volfil_vtkShape_180526.py] v0.1
    # OUT_FILE = 'cubeGroup_180526_t1040.stl'
    OUT_FILE = 'cubeGroup_180527_t0805.stl'
    DIPOLE_SIZE = 0.25

    with open(IN_FILE, 'rb') as rfp:
        dipoles = np.genfromtxt(IN_FILE)

    with open(OUT_FILE, 'wb') as wfp:
        writer = STLWR.Binary_STL_Writer(wfp)

        for adipole in dipoles:
            cube = get_cube(DIPOLE_SIZE, adipole)
            writer.add_faces(cube)

        writer.close()
    print('OUT:', OUT_FILE)

if __name__ == '__main__':
    make_cubeGroup()

結果

結果 20x20x20

20x20x20: 7分20秒の処理 (volfil_vtkShape_180526.py)。
331 KBのSTLファイル。

make_cubeGroup_180520.pyの処理は数秒。

cubeGroup_180527_t0805.stl @ GitHub

処理時間はかかるが、閲覧時は「XXソフトが必要」なく、ある程度の性能のあるPC上で閲覧できる。

結果 40x40x40

分解能を上げてみた。

40x40x40: 52分25秒の処理。
4894 dipoleとなった。

cubeGroup_180527_t0935.stl @ GitHub

生成したパラメータでのGaussin Random Sphere形状の特徴が捉えやすくなった。

結果 80x80x80

80x80x80: 380分24秒の処理。
40,662 dipoleとなった。

(注意: 下記のファイルは23.3MBです。)

cubeGroup_180527_t1052.stl @ GitHub