gnuplotの3Dグラフには少し癖がある．ここでは，3Dメッシュグラフを滑らかにする方法を記録する．
splot {数式}であればisosampleで間隔を調整できるが，数値計算の結果などsplot {datafile}だとisosampleが効かない．データファイルを間引いて描画すればいいのだが，そうすると曲線の解像度が落ちて美しくない．
そこで，多数の線を描画することにより面を構成する方法を記載する．
なお，めんどくさいので目盛りは省略した．
Microsoft Windows 10 Home (64bit)
Gnuplot Version 5.0 patchlevel 4
とりあえす，まずデータファイルを作成する．まあ，なんでもいいのだがここでは，
$$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 ...
というものである．なお，
<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データ構造を用いる．）ただし，この時の描画コマンドは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
で描画したのが下図である．
これだとメッシュが細かすぎてよくわからない．メッシュ数を減らしたかったのだが，isosampleなどが効かずにできなかった．
方法を知っている方はぜひコメントにて教えていただきたい．
章のデータファイルを間引いてメッシュを荒くした．データファイル変換スクリプトのソースコードは「」- , である．
set ticslevel 0 unset xtics unset ytics unset ztics splot "cnv1.dat" matrix with lines
で描画すると下図のようになる．特に最大値付近を見るとよく分かるが，格子点間を直線で描画しているため，なめらかな曲線とはいえない．
メッシュの線を曲線として描画させたのが，下図である．
set ticslevel 0 unset xtics unset ytics unset ztics splot "cnv2.dat" with lines
データ構造は，
<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>
である．ブロック間に２行空行をあけることと，x,yそれぞれに対して曲線を用意することがポイントである．
データファイル変換スクリプトのソースコードは「」- , である．
#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;
}
#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;
}
#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;
}
#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
名前
Email (※公開されることはありません)
コメント