geometry > Link > Shapely > lineとpolygonのintersection > trap:点でなく線分が返される場合がある

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

import matplotlib.pyplot as plt
from matplotlib import cm, colors
from mpl_toolkits.mplot3d import Axes3D
from pylab import rcParams
import numpy as np
import time
import datetime as dt
from shapely.geometry import LineString, Polygon

"""
v0.1 Apr. 30, 2018
  - cut with plane
=== branched from [showChebyshev_170910.ipynb] v0.4 ===
v0.4 Sep. 10, 2017
  - read ADDA file [IntField-Y]
v0.3 Sep. 10, 2017
  - set colors in X direction
  - increase number of spheres to 3000
v0.2 Sep. 10, 2017
  - show 8 spheres
  - lower the resolution of the sphere (from 100j to 6j)
v0.1 Sep. 10, 2017
  - show 2 spheres
"""

# coding rule: PEP8

rcParams['figure.figsize'] = 15, 10


# reference
# https://stackoverflow.com/questions/31768031/plotting-points-on-the-surface-of-a-sphere-in-pythons-matplotlib


def plot_spheres(xps, yps, zps):
    for elem in zip(xps, yps, zps):
        axp, ayp, azp = elem
        # print(elem)
        dx = x + axp
        dy = y + ayp
        dz = z + azp
        ax.plot_surface(
            dx, dy, dz,  rstride=1, cstride=1, color='c', 
            alpha=1.0, linewidth=0,
            facecolors=plt.cm.Set2((dx - 0) / (50 - 0)))  # 50: arbitrary chosen to set colors


start_time = time.time()

# Create a sphere
r = 1
pi = np.pi
cos = np.cos
sin = np.sin
phi, theta = np.mgrid[0.0:pi:6j, 0.0:2.0*pi:6j]
x = r*sin(phi)*cos(theta)
y = r*sin(phi)*sin(theta)
z = r*cos(phi)

RAD_OUTER = 10

# Cubdoid
NUM_RESOL = 20
inx = np.linspace(-11, 11, NUM_RESOL)
iny = np.linspace(-11, 11, NUM_RESOL)
inz = np.linspace(-11, 11, NUM_RESOL)
print(inx)
mx, my, mz = np.meshgrid(inx, iny, inz)

# Triangles to define the shape of the particle
DIST_POLY = 7
plys = [ Polygon([[0, DIST_POLY, 0], [DIST_POLY, 0, 0], [DIST_POLY, DIST_POLY, DIST_POLY]])]
# plys += [ Polygon([[0, DIST_POLY, DIST_POLY], [DIST_POLY, 0, DIST_POLY], [0, 0, DIST_POLY]])]

xpar = []
ypar = []
zpar = []
for idx in range(len(mx)):
    for idy in range(len(my)):
        for idz in range(len(mz)):
            inside = True
            for apoly in plys:
                #xs_ray, ys_ray, xs_ray = calc_raylinegeometry()
                a1 = np.array([0.0, 0.0, 0.0])
                a2 = np.array([mx[idx][idy][idz], my[idx][idy][idz], mz[idx][idy][idz]])
                sga = LineString([(a1[0], a1[1], a1[2]), (a2[0], a2[1], a2[2])])
                its = sga.intersection(apoly)
                if its:
                    inside = False

            if inside:
                xpar += [mx[idx][idy][idz]]
                ypar += [my[idx][idy][idz]]
                zpar += [mz[idx][idy][idz]]

# Set colours and render
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

SKIP_NUM = 1
xp = xpar[::SKIP_NUM]
yp = ypar[::SKIP_NUM]
zp = zpar[::SKIP_NUM]

print(dt.datetime.now())

plot_spheres(xp, yp, zp)

print(dt.datetime.now())

ax.set_xlim([-30, 30])
ax.set_ylim([-30, 30])
ax.set_zlim([-30, 30])
ax.set_aspect("equal")

ax.view_init(elev=30, azim=0)

plt.tight_layout()
plt.show()

duration = time.time() - start_time
print("%.3f" % duration)