日々適当

スクリプトでのカーブの描きかた？

import bpy
import math

coords = []

startAngle = 30.0
endAngle = 300.0
pointNo = 20

angleStep = (endAngle - startAngle) / ( pointNo - 1)

for i in range( 0 ,pointNo ):
    print( angleStep * i )
    coords.append( ( math.sin( ( angleStep * i  + startAngle )/ 360. * math.pi * 2. )  , 0.0,  math.cos( ( angleStep * i  + startAngle )/ 360. * math.pi * 2. )))
    
curveData = bpy.data.curves.new( 'myCurve', type='CURVE' )
curveData.dimensions = '3D'
curveData.resolution_u = 1

polyline = curveData.splines.new( 'POLY' )
polyline.points.add( len( coords ) - 1 )
for i, coord in enumerate( coords ):
    x, y, z = coord
    polyline.points[ i ].co = ( x, y, z, 1 )
    
curveOB = bpy.data.objects.new('myCurve', curveData)

scn = bpy.context.scene
scn.collection.objects.link(curveOB)
bpy.context.view_layer.objects.active = curveOB
curveOB.select_set(True)

270度の円弧を描いてみたものです。参考にしたサイト（ほぼそのまま）は Blender Stack Exchange です。
理解はしきってないけど、カーブを作って、それと関連づけたオブジェクトを作って、そのオブジェクトをシーンにリンクしたって流れらしいってことまでは分かった、ような気がする。

curveDataとpolyline、それぞれが作られた後、print( type( curveData ) ) や print( type( polyline ) ) ってしてみると

class 'bpy.types.Curve'
class 'bpy.types.Spline'

というお返事。
bpy.types.Splineクラスは

Element of a curve, either NURBS, Bezier or Polyline or a character with text objects

で、bpy.types.Curveは

Curve data-block storing curves, splines and NURBS

とドキュメントには書かれておりました。んー

Siggraph 2020が開催されていたのですね。知らんかった。で、そこで例年通りMAXONがイベントをネットで行っていたのが日本時間昨夜ということみたいです。
つまり、Redshift 3dのMac版が姿を見せました。MAXONのTwitterによると

昨晩発表されたRedshift関連のアナウンスです。
🎉Metal対応（数週間後パブリックベータ予定）
🎉Redshift RT
🎉Blenderのサポートhttps://t.co/KEXMibaGwz

— MaxonJapan (@maxonjapan) July 22, 2020

ってことみたいで、いよいよ使える日も近いって感じです。Metal版もBlenderで動作するようになるならますます嬉しいかも。です。

まぁファーストバージョンにそこまでは期待しないけど、これが順調にバージョンを重ねていけるよう、Appleにおかれましてもちゃんとサポートをして欲しいものです。

とかいう商業レベルのお話ではなく、世界は確実に進歩していると感じられるこーいうイベントは素晴らしいなぁと思うのは毎年のことでございます。

SIGGRAPH 2020: Technical Papers Preview Trailer

ほら、陰影までベイクしたテクスチャをそのまま表示するようなCGを作りたい時あるじゃないですか。そうするとベイクした結果を貼り付けてすぐに確認したくなるんですけど、Blenderだとそんな時どうするんでしょか。

こんな感じで陰影までベイクしてみたものの、それを貼り付けた状態を確認する時、シェーダをMix Shaderで分岐させて Fac パラメータを0か1かして、テクスチャのみの質感か元々のシェーダの質感かを切り替える。

そんな考え方でいいんでしょうか？

Boolean Logic的なノードはないみたいなんですけど、あるんすかね？
まぁともあれ、そんな感じで設定しておいて、どっかでその設定をしたマテリアルをためておいて一括で変更できるといいですよね。ってかそんな仕組みはないですかね。ベイクしたものを貼り付けた結果を確認したいときの王道ってどんなだろ。別に専用のシーンを持たせておいてそっちで確認するとか？

Blender 2.83のベータっす。

マテリアルとUV設定されたオブジェクトでカラー情報をテクスチャにベイクしようと苦闘しておりました（Bake TypeをDiffuseにし、Colorのみの影響された状態）。マテリアルはPrincipled BSDFを繋いだだけのシンプルな構造です。
何に苦闘していたかと言えば、ベイクの結果が

こんな感じで黒くなってしまうことに苦慮してまして、もともといじっていたものはこの赤い部分のマテリアルが黒色のものだったから真っ黒になって出てきていたわけです。ために、レンダリングの設定やオブジェクトの構造の問題かなと悶々としてたのだけど、マテリアルが原因の模様。

こんな具合のマテリアルなのだが、赤い部分と青い部分で何が違ったかというと、Metallicの値が違っておりました。赤いのはMetallicが0だったのです。
ってことで青い方もMetallicを0にしまして

ベイクしてみますと、

ちゃんと青色がのりました。
あーそっか、PBRだからMetallicをオンにしちゃうと、BaseColorはDiffuseには割り振られないってことか。これはうっかりでございます。

この地図に合わせてプロットすべしって指令がありました。プロット元のデータには緯度・経度が記されているから楽なもんだぜ、見た感じメルカトル図法だし、って思ったらミラー図法でした。

メルカトル図法の、南北両極が無限遠点になってしまうという問題を改善した図法で、地理緯度を 4/5 倍してからメルカトル図法で投影して、縦方向に 5/4 倍する。つまり地球を半径を1とする単位球とみなしたとき、ミラー図法において経度 λ, 地理緯度φから地図上の点 x, y へ投影する座標換算式は次式で与えられる：

ミラー図法 - Wikipedia

って書かれていました。説明の次式の中にgdの逆関数が混じっており「？」とね。グーデルマン関数の逆関数であると記されています。さらにWikipediaを漁るとグーデルマン関数の逆関数は最終的には arsinh( tan x ) なのだそうです。arsinh？ってなるんですよ。逆双曲線関数って言うんだそうです。arsinh z = ln( z + Sqrt( z^2 + 1 ) ) って書いてありました。あとは ln がどうにかなればどうにかなるね（自然対数が用意されていないスクリプト環境なもので）って思いつつ、なんか楽にできないものかなとさらに漁ると、lnは 2.303 * log x なのだそうです。この場合のlogは常用対数。常用対数は用意されているからこれでできるなってやってきたけど、なんかうまくいかない。<追記>と思ったらこの環境のLog()って自然対数のことだった。ってことでうまくいかない原因は2.303倍していたせいでした。</追記>
多分どこか間違っているんだろうけど時間が無いって状態なので、最初のWikipediaの文章を読み返すと「経度を4/5倍してメルカトル図法で投影して、縦方向に5/4する」って書いてあります。
これは単にメルカトル図法にできればいいだけじゃね？ってことでございました。メルカトル図法は sin の逆関数と tanh を導き出せればOKで、使用しようとしている環境では両方とも用意されているのでいけました。

ちなみにWikipediaを見るとメルカトル図法はグーデルマン関数だけで最終的には表現されておりました。最初のWikipediaの説明に納得した昨晩でございました。

先日のmatplotlibを使ってsvgに書き出す場合、日本地図全体の中のその県の位置に描画されてほしかった問題は、日本全体を書き出すならどんどん追記してやればいいだけなので、あまりmatplotlibをお勉強する必要はありませんでした。

これは行政区域データの1番から10番のZIpファイルから作ったものです。コードは以下。

import os
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import geopandas
from shapely.ops import cascaded_union

fig , ax = plt.subplots()

for i in range( 1 , 11 ):
	if i < 10:
		fileNo = '0' + str(i)
	else:
		fileNo = str(i)
	fileName = "N03-190101_" + fileNo + "_GML.zip"
	f = os.path.join( 'zip://pathToDirectory' , fileName)
	src = geopandas.read_file( f  , encoding="SHIFT-JIS", )
	#Web Mercator projection
	srcCnv = src.to_crs("EPSG:3857")
	polys = srcCnv.geometry
	polys = cascaded_union(polys)
	boundary = geopandas.GeoSeries( polys )
	boundary.plot( ax=ax,  color = cm.hsv( 1./47 * i  ))

plt.savefig('out_webMerc.svg')

とは言え、matplotlibのfigure, axes, axis なんてところをちゃんと理解せんといかんすね。

ということで、昨日のコードを変更してメルカトル図法にしたいと思いました。

GeoDataFrame.to_crs() ってコマンドを使うことで実現するようです。座標系の参考としては
List of map projections [Wikipedia]
とかみると良さそう。とりあえずWebメルカトルなEPSG:3857を指定してみました。

import os
import matplotlib.pyplot as plt
import geopandas
from shapely.ops import cascaded_union

f = 'pathTo/N03-19_190101.shp'
src = geopandas.read_file( f  , encoding="SHIFT-JIS", )
#Web Mercator projection
srcCnv = src.to_crs("EPSG:3857")
polys = srcCnv.geometry
polygons = [ poly for poly in polys ]
boundary = geopandas.GeoSeries(cascaded_union(polygons))

boundary.plot(color = 'red')

plt.savefig('out_webMerc.svg')

ってことでここまでは問題なくできたように思います。
で、問題はこれを各県とか政令指定都市とかで行ったとき、最終的にプロットされる範囲をどうしたもんか問題で、日本列島の一部とした範囲にプロットされて欲しいのだけど、上記コードだとそれができないのですね（そりゃそうなんだけど）。matplotlibのお勉強が必要です…

ちょっと必要に迫られてGeoPandasってのを試してる。インストールはこちら[GIS奮闘記] を参考にオプション含め、下のようにインストールを行いました。python3 -m venv venv を走らせたプロジェクトフォルダ内で作業してます。

[venv] [mac:pathToPrj] userName% pip3 install numpy
（略）
Installing collected packages: numpy
Successfully installed numpy-1.18.2
[venv] [mac:pathToPrj] userName% pip3 install pandas
（略）
Successfully installed pandas-1.0.3 python-dateutil-2.8.1 pytz-2019.3 six-1.14.0
[venv] [mac:pathToPrj] userName% pip3 install shapely
（略）
Successfully installed shapely-1.7.0
[venv] [mac:pathToPrj] userName% pip3 install fiona
（略）
Successfully installed attrs-19.3.0 click-7.1.1 click-plugins-1.1.1 cligj-0.5.0 fiona-1.8.13.post1 munch-2.5.0
[venv] [mac:pathToPrj] userName% pip3 install six
Requirement already satisfied: six in ./venv/lib/python3.7/site-packages (1.14.0)
[venv] [mac:pathToPrj] userName% pip3 install pyproj
（略）
Successfully installed pyproj-2.6.0
[venv] [mac:pathToPrj] userName% pip3 install descartes
（略）
Successfully installed cycler-0.10.0 descartes-1.1.0 kiwisolver-1.2.0 matplotlib-3.2.1 pyparsing-2.4.7
[venv] [mac:pathToPrj] userName% pip3 install matplotlib
Requirement already satisfied: matplotlib in ./venv/lib/python3.7/site-packages (3.2.1)
（略）
Requirement already satisfied: six in ./venv/lib/python3.7/site-packages (from cycler>=0.10->matplotlib) (1.14.0)
[venv] [mac:pathToPrj] userName% pip3 install geopy
（略）
Successfully installed geographiclib-1.50 geopy-1.21.0
[venv] [mac:pathToPrj] userName% pip3 install geopandas
（略）
Successfully installed geopandas-0.7.0
[venv] [mac:pathToPrj] userName% 

実際に書いてみたところ、Spatialindexってのが足りなかったんで、これはbrewで入れてやりました。
brew install Spatialindex
空間インデックスって概念のものを取り扱うライブラリってことでしょか。

これで国土地理院国土数値情報の行政区域データから日本列島とか千葉県とかそんな塊の形状を出したかったわけです。
行政区域データは市区町村の形状が入っているために、千葉県の形が欲しい場合とかに市区町村境界部分が無駄なのです。じゃぁIllustratorでマージすればいいじゃんって思うのですけど、Illustratorは一つの形状が32000ポイントを超えることができないため、形状をマージしていく過程でそれの数を超えてしまう事が多々ありどうにもならない。それ以上のポイントを取り扱えるAfinity Designerだと、きれいにマージしてくれない（境界線部分にゴミが残ることが多い）のです。だからPythonで処理できないかなと、形状をマージすることができる何かを探して、GeoPandasに行き当たりました。

というわけで、各所参考に組んでみたのがこちら。（合成部分はこちらを参考にしました。）

import os
import matplotlib.pyplot as plt
import geopandas
from shapely.ops import cascaded_union

f = 'pathTo/N03-190101_GML/N03-19_190101.shp'
#行政地域データの日本全土のもの
src = geopandas.read_file( f  , encoding="SHIFT-JIS", )
polys = src.geometry
polygons = [ poly for poly in polys ]
boundary = geopandas.GeoSeries(cascaded_union(polygons))

boundary.plot(color = 'red')
plt.savefig('out.svg')

みた感じ、割といい感じの結果が出たような気がします。

上は全体、下は犬吠埼あたりを拡大したもの。
これ、図法変換をうまく行うことができれば、以前やってた市区町村区切りを持った形の日本地図を出す場合でも、以前よりだいぶシンプルに処理できるかもしれないと期待が高まってきております。

ということで、形状を維持した形で点の数を間引きたくなりました。
38000個の点で構成された線を32000個の点にまで減少させたい時、6000個減らしたいので、38000/6000=6.33...ってことで7個おきに間引けば、一応32000個を下回らせることが可能ですけど、当然形状はよろしくないことになります（事実上あまり気にならないかもしれないけどね）。

ということで、形状をできるだけ維持した形での点を間引く方法はなんだろうと調べると、Douglas and Peuckerアルゴリズムってのが有名なのだそうです。点と点を直線で結んで、それらの間の点とその直線（傾向線）の距離が許容値（ε）より小さい場合にその点を削除し、また次の直線を引いて点との距離が許容値より小さければ点を削除する、というのを削除する点がなくなるまで繰り返す手法は、形状を特徴づける重要な点が維持されるという利点がある一方、見栄えがちょと悪化することがあるようです。まぁでも、単純に機械的に点を削減するよりはマシだろうと考え、この方法を採用する時、Pythonのモジュールとして rdp [GitHub]ってのが配布されていました。普通に pip3 install rdp で入ります。

ちなみに、元の線の曲率を維持するようにして点の数を減らす手法もあり、こっちは形はいい感じになるけど点の位置はオリジナルから外れますね。それをできるモジュールは知りません。OpenCVとか使うのかね？

rdpの基本的な使い方は単純でした。座標の入った配列と許容値 epsilon を指定するだけ。許容値は傾向線からの最小距離ですね。だから値が小さいほど消される点の数は少なくなります。

import svgwrite
import xml.etree.ElementTree as ET
import numpy as np
from rdp import rdp

#二つの点の距離を返す
def distance( pos0 , pos1 ):
	p0 = np.array( pos0 )
	p1 = np.array( pos1 )
	v = p1 - p0
	return np.linalg.norm( v )

#元となる形状の入ったsvgファイルを読み込みその座標を配列 poses に入れる。
tree = ET.parse( 'testShape.svg' )
root = tree.getroot()
poses = []
for child in root:
	if 'path' in child.tag:
		posInfo = child.attrib[ 'd' ]
		for posStr in posInfo[1:-1].split( 'L' ):
			pos = posStr.split( ',' )
			poses.append( [ float( pos[0] ),   float( pos[1] ) ])

#ポイント数を80%に削減することを目標とすることにして、目標のポイント数を設定
newPointsNum = int( len( poses ) * 0.8)
baseLen = distance(  poses[ 0 ], poses[ 1 ] )

#目標の点数を下回るまで繰り返す。
#ここで問題の rdp モジュールを使ってる
i = 1
while True:
	newPoses = []
	eps = baseLen * i * 0.005 #距離の許容値の指定
	newPoses = rdp( poses , epsilon = eps )
	print( len( newPoses ) )
	if len( newPoses ) < newPointsNum:
		print( '結果 ', len( newPoses ) )
		break
	i += 1

#svgファイルのキャンパスサイズを決めるためのとこ
pMaxX = max( poses, key = lambda x: x[ 0 ]  )
pMinX = min( poses, key = lambda x: x[ 0 ]  )
pMaxY = max( poses, key = lambda x: x[ 1 ]  )
pMinY = min( poses, key = lambda x: x[ 1 ]  )
w = pMaxX[ 0 ] - pMinX[ 0 ] + pMinX[ 0 ] * 2
h =  pMaxY[ 1 ] -  pMinY[ 1 ]  + pMinY[ 1 ] * 2

#svgファイルへの書き出し。
dwg = svgwrite.Drawing( 'result_testShape.svg' , size=( w , h))
dwg.add( dwg.polygon( points=poses ,stroke='blue', stroke_width=8.0, fill='none', id='オリジナル'))
dwg.add( dwg.polygon( points=newPoses ,stroke='red', stroke_width=4.0, fill='none', id='結果' ) )
dwg.save()

得られた結果が以下。青線がオリジナルの形状。赤線が点の数を削減したもの。

目標とする点の数を指定することができないので、とりあえず目標の点の数になるまでループを回しているのだけど、点数が多いと時間がかかりそうでどうしたもんか…。あとepsilonの値を適切に与えるためにもどうしたらいいのか悩みどころっす。

SVGのパスデータにおいて、一つのパスオブジェクトは何個の点まで持てるかの確認です。ググると32000ってことらしいのですが、その検証です。

import svgwrite
import math

pointNum = 100
p = ""
for i in range( pointNum ):
	x = 100.0  * math.cos( math.radians( 360.0/float( pointNum ) * i ) ) + 105.0
	y = 100.0  * math.sin( math.radians( 360.0/float( pointNum ) * i ) ) + 105.0
	if i == 0:
		p += " M " + str( x ) + " " + str( y )
	elif i == pointNum -1:
		p += " L " + str( x ) + " " + str( y ) + "Z"
	else:
		p += " L " + str( x ) + " " + str( y ) 
dwg = svgwrite.drawing.Drawing(  'test_path.svg'  , size=( 210, 210), profile='full' )
path = dwg.path( p ).stroke( color="red" , width=0.5 ).fill( "none" )
dwg.add( path )
dwg.save()

こんな感じのスクリプトで指定したポイント数の円を描いたsvgファイルをIllustratorに読み込みます。この場合、pointNumの値が100なので100ポイントで構成された円が出力されます。
このpointNumの値を32000にしたものはちゃんとIllustratorに読み込めました。しかして、これを32001にするとダメです。上限は32000ということみたいです。

ちなみにですが、先のスクリプトのfor文部分をもう一個追加して、その中の円の半径の指定部分を小さくし、つまり二重の円を描いた時、pointNumが32000でも読み込まれます。ポイント数自体は64000になりますが、一つ一つの円は32000なのでOKってことですね。。

svg経由でこの制限を突破するのは難しいですかねー？

<追記>そもそも、Illustratorの仕様としてパス一個につき32000ポイントしか許容していない疑惑。であるならポイント数を削減する作業が必要になるのだろうけど、隣接する区域のパスとズレが生じそうで嫌やんだよなぁ…。ちなみに先のエントリで話題にした佐渡市の佐渡島本島部分のポイント数は37662でございました。</追記>