2 次元データの補間

2010年1月13日

はじめに

ParaView を使って、VTK の 2 次元点データ列から、CSV 形式の任意の 2 次元点列に値を補間する。

バージョン

ParaView 3.6.2

ファイル

• contour.vtk.zip
• point.csv

サンプルとしてお使いください。VTK ファイルのほうは ZIP 圧縮しているので、展開してください。

任意の点におけるデータの取得

ParaView でデータを表示させたとき、そこから任意の点でのデータを得るには、フィルタの Probe Location を使う。

メニューの Filters から Probe Location を選ぶと、点の座標を入れるところがあるので、そこに座標値を入れて Apply ボタンを押すと、その座標値におけるデータが得られる。マウスポインタを画面上に置いてキーボードの 'p' を押すことで、ポインタが指している座標値を設定することもできる。

具体的な値を知るには、Information タブの Data Arrays を見る。レンジで表示されてあるが、1 点分のデータしかないので、最小値も最大値も同じである。あるいは、画面右上の小さいボタンで画面を分割して、Spreadsheet View で見るという手もある。

ちなみに、Number of Points と Radius を設定すると、指定した座標周りに設定した半径で点を適当に配置してくれる。

この機能を使えば、2 次元点データをデローニー分割した後、任意の点での値を取得することができる。とはいえ、値を取りたい位置がたくさんあると、いちいち手でやってはいられない。Python での自動化を試みよう。

Python による自動化

Python スクリプトを示す。

vtkToPoint.py

# vtkToPoint

import os

def vtkToPoint(vtkFilename, csvFilename):
	reader = LegacyVTKReader(FileNames=vtkFilename)
	view = GetActiveView()
	view.StillRender()
	
	d = Delaunay2D(reader)
	view.StillRender()
	
	p = ProbeLocation()
	Show(p)
	
	pt = p.ProbeType
	data = p.PointData['point_scalars']
	
	inFile = open(csvFilename, "r")
	inFile.readline()
	
	vtkBasename = os.path.splitext(vtkFilename)[0]
	csvBasename = os.path.splitext(csvFilename)[0]
	
	outFile = open(vtkBasename + '-' + csvBasename + '.csv', 'w')
	outFile.write('x,y,value\n')
	
	for line0 in inFile:
		line = line0[0:len(line0)-1]
		item = line.split(",")
		x = float(item[0])
		y = float(item[1])
		
		pt.Center = (x, y, 0.)
		view.StillRender()
		
		v = data.GetRange()[0]
		
		outFile.write('%f,%f,%f\n' % (x, y, v))
	
	inFile.close()
	outFile.close()

使い方は、まずメニュー Tools から Python Shell を開き、Run Script ボタンで上のスクリプトを読み込む。そして、以下のように関数を呼び出す。

>>> vtkToPoint('contour.vtk', 'point.csv')

第 1 引数がデータの VTK ファイル、第 2 引数が値を取りたい点列が書かれた CSV ファイルである。上の場合、'contour-point.csv' という点データ列の CSV ファイルが書き出される。入力する点列の CSV ファイルの書式は、以下のように 1 行目が項目のタイトル、それに続いて 2 次元座標値が並んでいる。

x,y
1.374302387,3.228690147
-4.621094227,-11.71187782
-2.019116879,2.14136982
...

作業をするにあたって、作業フォルダをカレントのフォルダにするために、Windows 版であれば、ParaView のアイコンを作業フォルダにコピーし、アイコンを右クリック、プロパティでショートカットタブの作業フォルダを空欄にし、そのアイコンで ParaView を起動するとよい。ただし、ヘルプを開けなくなるが。

スクリプトの内容を説明する。

	reader = LegacyVTKReader(FileNames=vtkFilename)
	view = GetActiveView()
	view.StillRender()

LegacyVTKReader() で VTK ファイルを読み込む。現在のビューを取得して、StillRender() を呼び出す。これが何なのか知らないが、これがないとデータが正しく取れない。以下、ことあるごとに呼び出す。

	d = Delaunay2D(reader)
	view.StillRender()

フィルタの Delaunay2D の適用。

	p = ProbeLocation()
	Show(p)

フィルタの Probe Location を適用する。ここで表示する。データを取るだけなので、わざわざ表示しなくてもよさそうなものだが、こうしないとちゃんとデータが取れない。

	pt = p.ProbeType
	data = p.PointData['point_scalars']

後で使うために ProbeLocation の ProbeType の参照を得る。また、取り出すデータとして、'point_scalars' を指定している。もし別のデータを取得したければ、このあたりを書き換えること。

その後、入力用と出力用のファイルを開いて、座標値を読み込んで、以下でデータを取る座標を指定している。

		pt.Center = (x, y, 0.)
		view.StillRender()

そして、値の取得。

		v = data.GetRange()[0]

レンジから値を得ている。最小値と最大値は同じ値なので、最小値を取ってその点の値としている。

入力データと出力データを比べたのが以下の絵である。左が元のデータ、右が補間されたデータ。

点の数が少ないのでちょっとぼんやりして見えるが、まあよさそうである。