0
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

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

Last updated at Posted at 2018-04-29
動作環境
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(星形)のedgeの交点を見つけて、線分を引く (平面角すべての範囲で行う) > find_rayCrossing_starShapedEged_180417 v0.5

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

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

処理方法

下記としてみた。

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

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

code v0.1

findDipolesInisideStarShape_180429.ipynb
%matplotlib inline

import numpy as np
import sympy as sp
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

'''
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 = 20
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())

for idx in range(len(mx)):
    print(idx, end='.')
    for idy in range(len(my)):
        # segment from center to dipole
        pa1 = sp.Point(0.0, 0.0)
        pa2 = sp.Point(mx[idx][idy], my[idx][idy])
        sga = sp.Segment(pa1, pa2)
        # edge of the shape
        isinside = True
        for aidx in combs:
            if not isinside:
                continue
            if list(aidx) not in edgeidx:
                continue
            # b1 = np.array([xs_str[aidx[0]], ys_str[aidx[0]]])
            # b2 = np.array([xs_str[aidx[1]], ys_str[aidx[1]]])
            sgb = get_linesegment(xs_str, ys_str, aidx[0], aidx[1])
            its = sp.intersection(sga, sgb)
            if its:  # outside
                isinside = False
                continue
        if isinside:
            ax1.scatter(mx[idx][idy], my[idx][idy])
   

print(dt.datetime.now())

fig.tight_layout()

qiita.png

処理に3分かかる。

3次元にした場合、膨大な待ち時間になる。

0
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?