動作環境

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

前回: geometry > star shape(星形)の内部をdipoleでfillする > findDipolesInisideStarShape_180429 v0.1 > sympy使用 > 3分かかる (MESH_RESOL=20)

Gaussian random sphereの形状情報をもとにADDAで光散乱数値シミュレーションを行うには、dipoleでfillした形状に変換する必要がある。

星形をもとに方法を検討してきた。

処理方法

下記としてみた。

格子状のdipoleを用意する
中心からdipoleの線分を検討する。
線分が形状の線分と交差しているか
- 交差あり: 外側
- 交差なし: 内側

内側のdipoleだけを残すことで、形状をdipoleでfillした状態になる。

code v0.2

前回の処理(sympy使用)ではMESH_RESOL=20に対して３分かかる。

遅い。

sympyのintersectionチェックは自分が使う機能以上に処理をしているように思われる (線分同士の交差だけでなく、線分と多角形の交差など)。
もっと軽量な処理を探してみた。

shapelyを使った実装例を以下に見つけた。
https://rosettacode.org/wiki/Find_the_intersection_of_two_lines#Python

使ってみた。

findDipolesInisideStarShape_180429.ipynb

%matplotlib inline

import numpy as np
import sys
import matplotlib.pyplot as plt
from pylab import rcParams
from itertools import combinations
import geometry_starShaped_180428 as GSS
import datetime as dt
from shapely.geometry import LineString

'''
v0.2 Apr. 29, 2018
  - use shapely instead of sympy
v0.1 Apr. 29, 2018
  - use algorithm to check inside/outside the given shape
      + segment (from center to dipole and segment of the shape
      + crossing segments indicates dipole is located outside the shape
=== branched from [find_rayCrossing_starShapedEged_180417.ipynb] ===
v0.5 Apr. 29, 2018
  - loop through [theta_deg]
  - add get_linesegment()
  - add calc_raylinegeometry()
v0.4 Apr. 29, 2018
  - find crossing point instead of crossing edge
     + remove: import [geometry_lineintersect_180415]
     + import sympy
v0.3 Apr. 28, 2018
  - import geometry_starShaped_180428 (v0.3)
v0.2 Apr. 28, 2018
  - exclude vertices pairs not included in the output from [geometry_starShaped_180428]
  - import geometry_starShaped_180428 (v0.2)
  - remove: import geometry_starShaped_180415
v0.1 Apr. 17, 2018
  - check the crossing   
      + import [geometry_lineintersect_180415]
  - define the ray from the center
    + draw the ray
    + add [xs_ray], [ys_ray]
    + add [theta_deg]
  - add idxs_seq[], combs[]
  - import [combinations]
  - branched from [geometry_starShaped_180414.ipynb]
'''

rcParams['figure.figsize'] = 14, 7
rcParams['figure.dpi'] = 110


def calc_raylinegeometry(theta_deg):
    # define the ray from the center
    radius = RAD_OUTER * 1.1  # 1.1 to cross the outmost edge
    xs_ray, ys_ray = [0.0], [0.0]  # center
    xs_ray += [radius * np.cos(np.deg2rad(theta_deg))]  # outmost
    ys_ray += [radius * np.sin(np.deg2rad(theta_deg))]  # outmost
    return xs_ray, ys_ray


def get_linesegment(xs, ys, idx0, idx1):
    pa1 = sp.Point(xs[idx0], ys[idx0])
    pa2 = sp.Point(xs[idx1], ys[idx1])
    return sp.Segment(pa1, pa2)

MESH_RESOL = 50
inx = np.linspace(-11, 11, MESH_RESOL, endpoint=True)
iny = np.linspace(-11, 11, MESH_RESOL, endpoint=True)
mx, my = np.meshgrid(inx, iny)

# 1. obtain star shaped points
RAD_INNER = 5
RAD_OUTER = 10
xs_str, ys_str, edgeidx = GSS.get_starShaped(RAD_INNER, RAD_OUTER)

# 2. obtain combinations of points
idxs_seq = range(len(xs_str))  # sequential indices to obtain combinations
combs = []  # combinations to obtain all the edges
for acomb in combinations(idxs_seq, 2):
    combs += [acomb]
print(combs)


fig = plt.figure()

ax1 = fig.add_subplot(1, 2, 1)
ax1.scatter(xs_str, ys_str, marker='*')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.grid(True)

print(dt.datetime.now())

cnt = 0  # counter to know the number of intersection check
for idx in range(len(mx)):
    print(idx, end='.')
    for idy in range(len(my)):
        # segment from center to dipole
        a1 = np.array([0.0, 0.0])
        a2 = np.array([mx[idx][idy], my[idx][idy]])
        sga = LineString([(a1[0], a1[1]), (a2[0], a2[1])])
        # edge of the shape
        isinside = True
        for aidx in combs:
            if not isinside:
                continue
            if list(aidx) not in edgeidx:
                continue
            cnt += 1
            b1 = np.array([xs_str[aidx[0]], ys_str[aidx[0]]])
            b2 = np.array([xs_str[aidx[1]], ys_str[aidx[1]]])
            sgb = LineString([(b1[0], b1[1]), (b2[0], b2[1])])
            if sga.intersection(sgb):
                isinside = False
                continue
        if isinside:
            ax1.scatter(mx[idx][idy], my[idx][idy])
   

print(cnt)
print(dt.datetime.now())

fig.tight_layout()

MESH_RESOL = 50でも１秒で処理できるようになった。

TODO

calc_raylinegeometry()使っていない。削除してもいい。