MENU

溶けかけてるうさぎ HP BLOG TOP RECENT ARTICLES POPULAR ARTICLES ABOUT THIS BLOG

CATEGORY

大学 (85) 航空宇宙 (55) 写真 (25) 旅行 (14) 飯・酒 (11) コンピュータ (88) その他 (13)

TAG

ARCHIVE

2018 (92) 2017 (80) 2016 (0)

RECENT

【駅メモ】4年目に突入して,ようやく3000駅突破 【WebRTC】Raspberry Pi搭載ロボットをWebRTCで遠隔操作しようとして失敗した 【航空宇宙】航空宇宙アドベントカレンダー 始まります! 【Perl】YAPC::Tokyo 2019 のチケットを確保しました! 【カメラ】Canonから富士フイルムに乗り換えました

【gnuplot】数値計算結果の3Dメッシュグラフを滑らかにする方法

2017-03-28

gnuplotの3Dグラフには少し癖がある.ここでは,3Dメッシュグラフを滑らかにする方法を記録する.

splot {数式}であればisosampleで間隔を調整できるが,数値計算の結果などsplot {datafile}だとisosampleが効かない.データファイルを間引いて描画すればいいのだが,そうすると曲線の解像度が落ちて美しくない.

そこで,多数の線を描画することにより面を構成する方法を記載する.

なお,めんどくさいので目盛りは省略した.

1.実行環境

Microsoft Windows 10 Home (64bit)

Gnuplot Version 5.0 patchlevel 4

.一般的なmatrixデータ構造の場合

とりあえす,まずデータファイルを作成する.まあ,なんでもいいのだがここでは,

$$z=\frac{\sin{\sqrt{x^2+y^2}}}{\sqrt{x^2+y^2}}$$

を描画しよう.

 

matrixデータ構造とは,

z11 z12 z13 z14 ...
z21 z22 z23 z24 ...
z31 z32 z33 z34 ...
matrixデータ構造

というものである.なお,

<N+1>  <y0>   <y1>   <y2>   ... <yN>
<x0>   <z0,0> <z0,1> <z0,2> ... <z0,N>
<x1>   <z1,0> <z1,1> <z1,2> ... <z1,N>
  :      :      :      :    ...   :
非一様matrixデータ構造

という非一様matrixデータ構造を用いれば,後から目盛り情報を付加する必要がなくて便利である.(ここでは一般的なmatrixデータ構造を用いる.)ただし,この時の描画コマンドはsplot {datafile} nonuniform matrix with linesとなる.

 

」に示す, で,データファイル(data_raw.dat)を作成して,

set ticslevel 0
unset xtics
unset ytics
unset ztics
splot "data_raw.dat" matrix with lines
gnuplot コマンド

で描画したのが下図である.

gnuplot splot matrix

これだとメッシュが細かすぎてよくわからない.メッシュ数を減らしたかったのだが,isosampleなどが効かずにできなかった.

方法を知っている方はぜひコメントにて教えていただきたい.

.matrixデータ構造を間引いた場合

章のデータファイルを間引いてメッシュを荒くした.データファイル変換スクリプトのソースコードは「」- , である.

set ticslevel 0
unset xtics
unset ytics
unset ztics
splot "cnv1.dat" matrix with lines
gnuplot コマンド

で描画すると下図のようになる.特に最大値付近を見るとよく分かるが,格子点間を直線で描画しているため,なめらかな曲線とはいえない.

gnuplot splot matrix

4.面描画ではなく,多数の線で面を構成する方法

メッシュの線を曲線として描画させたのが,下図である.

set ticslevel 0
unset xtics
unset ytics
unset ztics
splot "cnv2.dat" with lines
gnuplot コマンド
gnuplot splot

章の図に比べて美しくなっている.

データ構造は,

<x0>  <y0>  <z00>
<x0>  <y1>  <z01>
<x0>  <y2>  <z02>
  :     :     :
<x0>  <yN>  <z0N>


<x1>  <y0>  <z10>
<x1>  <y1>  <z11>
<x1>  <y2>  <z12>
  :     :     :
<x1>  <yN>  <z1N>

...

<xN>  <y0>  <zN0>
<xN>  <y1>  <zN1>
<xN>  <y2>  <zN2>
  :     :     :
<xN>  <yN>  <zNN>


<x0>  <y0>  <z00>
<x1>  <y0>  <z10>
<x2>  <y0>  <z20>
  :     :     :
<xN>  <y0>  <zN0>


<x0>  <y1>  <z01>
<x1>  <y1>  <z11>
<x2>  <y1>  <z21>
  :     :     :
<xN>  <y1>  <zN1>

...

<x0>  <yN>  <z0N>
<x1>  <yN>  <z1N>
<x2>  <yN>  <z2N>
  :     :     :
<xN>  <yN>  <zNN>
ここで使用するデータ構造

である.ブロック間に2行空行をあけることと,x,yそれぞれに対して曲線を用意することがポイントである.

データファイル変換スクリプトのソースコードは「」- , である.

5.ソースコード

#include "MyUtil.h"

using namespace std;
using namespace MyUtil;
using MyUtil::status_t;

const string OUTPUT_FILE_NAME = "data_raw.dat";
const double RANGE = 10.0;
const double dx = 0.1;
const int NX = static_cast<int>(2*RANGE / dx) +1;

double fn(double x, double y) {
	if (x==0.0 && y==0.0) return 1.0;
	return sin(sqrt(x*x+y*y)) / sqrt(x*x+y*y);
}

// 配列idxを座標に変換
double idx2x(double idx) {
	int centerIdx = NX / 2;
	return (idx - centerIdx) * dx;
}

int main() {
	ios::sync_with_stdio(false);
	vector<vector<double>> z(NX, vector<double>(NX, 0));
	for(int i=0; i<NX; ++i) {
		for(int j=0; j<NX; ++j) {
			z[i][j] = fn(idx2x(i),idx2x(j));
		}
	}

	output2dVector2Csv(OUTPUT_FILE_NAME, z, "\t");
	return 0;
}
Source.1 main.cpp(データファイル生成)
#include "MyUtil.h"

using namespace std;
using namespace MyUtil;
using MyUtil::status_t;

const string OUTPUT_FILE_NAME = "data_cnv1.dat";
const string INPUT_FILE_NAME = "data_raw.dat";
const int NESH_NUM = 21;

int main() {
	ios::sync_with_stdio(false);
	status_t status;
	vector<vector<double>> z;

	tie(status, z) = inputCsv22dVector<double>(INPUT_FILE_NAME, "\t");

	const int sizeX = z.size();
	const int sizeY = z[0].size();
	const int stepX = (sizeX-1) / (NESH_NUM-1);
	const int stepY = (sizeY-1) / (NESH_NUM-1);

	// z を間引いて出力.
	ostringstream oss;
	for (int i=0; i<sizeX; i+=stepX) {
		for (int j=0; j<(sizeY-1); j+=stepY) {
			oss << z[i][j] << "\t";
		}
		oss << z[i][sizeY-1] << "\n";
	}

	outputFile(OUTPUT_FILE_NAME, oss);
	return 0;
}
Source.2 convert_1.cpp(matrixデータ構造を間引くスクリプト)
#include "MyUtil.h"

using namespace std;
using namespace MyUtil;
using MyUtil::status_t;

const string OUTPUT_FILE_NAME = "data_cnv2.dat";
const string INPUT_FILE_NAME = "data_raw.dat";
const int NESH_NUM = 21;

int main() {
	ios::sync_with_stdio(false);
	status_t status;
	vector<vector<double>> z;

	tie(status, z) = inputCsv22dVector<double>(INPUT_FILE_NAME, "\t");

	const int sizeX = z.size();
	const int sizeY = z[0].size();
	const int stepX = (sizeX-1) / (NESH_NUM-1);
	const int stepY = (sizeY-1) / (NESH_NUM-1);

	// 目盛り用(今回は不使用)
	double dx = 0.1;
	double dy = 0.1;
	// 縦と横の線を独立して,面ではなく複数の線として出力
	outputFile(OUTPUT_FILE_NAME, Matrix2Gnuplot3DGridData(z, dx, dy, stepX, stepY));
	return 0;
}
Source.3 convert_2.cpp(多数の線によって面を構成するデータ構造へ変換するスクリプト)
#ifndef MYUTIL_H
#define MYUTIL_H

#define _USE_MATH_DEFINES
#include <cmath>
#include <iostream>
#include <sstream>      // ss
#include <fstream>      // ifstream, ofstream
#include <iomanip>      // マニピュレータ
#include <string>
#include <vector>
#include <tuple>
#include <time.h>

namespace MyUtil {

using namespace std;
using status_t = int;       // 成功失敗などの返り値用
const status_t STATUS_SUCCESS = 1;
const status_t STATUS_FAILURE = 0;

// 汎用Print
template <typename T> void p(T a);
template <typename T> void pn(T a);
// 文字列処理系
vector<string> SplitString(const string& input, const string delimiter);
template <typename T> vector<T> CastStringVector(const vector<string>& input);
// Util
tuple<status_t, vector<string>> inputFile2LineVector(const string& filename);
template <typename T> tuple<status_t, vector<vector<T>>> inputCsv22dVector(const string& filename, const string delimiter = ",");
inline status_t outputFile(const string& filename, const string& str, int isApp = 0);
inline status_t outputFile(const string& filename, const ostringstream& oss, int isApp = 0);
template <typename T> status_t output2dVector2Csv(const string& filename, vector<vector<T>>& vv, const string delimiter = ",");
template <typename T> string Matrix2Gnuplot3DGridData(vector<vector<T>>& vv, double dx, double dy, int stepX, int stepY, double offsetX=0, double offsetY=0);
void _PrintErr(const string& s);


// ########################################
// 汎用Print
template <typename T>
void p(T a) {
	cout << a << flush;
}

template <typename T>
void pn(T a) {
	cout << a << endl;
}

// ########################################
// 文字列処理系
vector<string> SplitString(const string& input, const string delimiter) {
	vector<string> result;
	if (input.size() == 0) return result;
	// delimiter = "" の場合
	if (delimiter == "" || delimiter.size() == 0) {
		for (int i=0;i<input.size();++i) {
			result.push_back(input.substr(i,1));
		}
		return result;
	}
	string::size_type pos = 0;
	while(pos != string::npos) {
		string::size_type p = input.find(delimiter, pos);
		if(p == string::npos) {
			result.push_back(input.substr(pos));
			break;
		} else {
			result.push_back(input.substr(pos, p - pos));
		}
		pos = p + delimiter.size();
	}
	return result;
}

template <typename T>
vector<T> CastStringVector(const vector<string>& input) {
	int size = input.size();
	vector<T> result(size);
	for (int i=0;i<size;++i) {
		istringstream iss(input[i]);
		iss >> result[i];
	}
	return result;
}

// ########################################
// Util
tuple<status_t, vector<string>> inputFile2LineVector(const string& filename) {
	vector<string> result(0);
	ifstream ifs;
	ifs.open(filename);
	string inputString;
	if (ifs.fail()) {
		_PrintErr("opne file error at inputFile2LineVector()");
		return forward_as_tuple(STATUS_FAILURE, result);
	}
	while (getline(ifs, inputString)) {
		result.push_back(inputString);
	}
	ifs.close();
	return forward_as_tuple(STATUS_SUCCESS, result);
}

template <typename T>
tuple<status_t, vector<vector<T>>> inputCsv22dVector(const string& filename, const string delimiter) {
	vector<vector<T>> result(0);
	status_t status;
	vector<string> stringLines(0);
	tie(status, stringLines) = inputFile2LineVector(filename);
	if (status == STATUS_FAILURE) {
		return forward_as_tuple(STATUS_FAILURE, result);
	}
	for(auto stringLine : stringLines) {
		result.push_back(CastStringVector<T>(SplitString(stringLine, delimiter)));
	}
	return forward_as_tuple(STATUS_SUCCESS, result);
}

inline status_t outputFile(const string& filename, const string& str, int isApp) {
	ofstream ofs;
	if (isApp == 1) {
		ofs.open(filename, ios::app);
	} else {
		ofs.open(filename);
	}
	if(! ofs.is_open()) {
		return STATUS_FAILURE;
	}
	ofs << str;
	ofs.close();
	return STATUS_SUCCESS;
}

inline status_t outputFile(const string& filename, const ostringstream& oss, int isApp) {
	return outputFile(filename, oss.str(), isApp);
}

template <typename T>
status_t output2dVector2Csv(const string& filename, vector<vector<T>>& vv, const string delimiter) {
	ostringstream oss;
	oss.precision(14);
	int rows = vv.size();
	int cols = vv[0].size();
	for (int i=0; i<rows; ++i) {
		for (int j=0; j<cols-1; ++j) {
			oss << vv[i][j] << delimiter;
		}
		oss << vv[i][cols-1] << "\n";
	}
	return outputFile(filename, oss);
}

template <typename T>
string Matrix2Gnuplot3DGridData(vector<vector<T>>& vv, double dx, double dy, int stepX, int stepY, double offsetX, double offsetY) {
	stringstream ss;
	int NX = vv.size();
	int NY = vv[0].size();
	for (int i=0; i<NX; i+=stepX) {
		for (int j=0; j<NY; ++j) {
			ss << i*dx+offsetX << "\t" << j*dy+offsetY << "\t" << vv[i][j] << "\n";
		}
		ss << "\n\n";
	}
	for (int i=0; i<NY; i+=stepY) {
		for (int j=0; j<NX; ++j) {
			ss << j*dx+offsetX << "\t" << i*dy+offsetY << "\t" << vv[j][i] << "\n";
		}
		ss << "\n\n";
	}
	return ss.str();
}

void _PrintErr(const string& s) {
	p("Error : ");
	pn(s);
	pn("Please press Enter key to continue...");
	getchar();
}

}
#endif
Source.4 MyUtil.h

6.関連記事

コメントを投稿

名前

Email (※公開されることはありません)

コメント