Instantly share code, notes, and snippets.

# komasaru/calc.cpp

Last active November 10, 2022 11:22
Show Gist options
• Save komasaru/fb31c087c71292b3a815535b8aae0fa1 to your computer and use it in GitHub Desktop.
C++ source code to interpolate by 3D-Spline.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 #include "calc.hpp" #include #include #include #include #include /** * @brief 配列サイズ取得 * * @return 配列サイズ(unsigned int) */ unsigned int Calc::get_size() { return data_src.size(); } /** * @brief 3次スプライン補間 * * @param[ref] n x 2 行列（配列） data_src (double) * @param[ref] ? x 2 行列（配列） data_res (double) * @return 真偽(bool) * @retval true 成功 * @retval false 失敗 */ bool Calc::interpolate_spline( std::vector>& data_src, std::vector>& data_res ) { std::vector a; // 配列 a std::vector b; // 配列 b std::vector c; // 配列 c std::vector d; // 配列 a std::vector h; // 配列 h std::vector v; // 配列 v std::vector w; // 配列 w std::vector> mtx; // 計算用行列 std::size_t i; // loop インデックス double s; // step double x; // 計算対象 x 値 double y; // 計算結果 y 値 unsigned int idx; // 計算対象インデックス try { // 配列 h 計算 if (!calc_h(data_src, h)) { std::cout << "[ERROR] H array could not be calculated!" << std::endl; return false; }; // 配列 w 計算 if (!calc_w(data_src, h, w)) { std::cout << "[ERROR] W array could not be calculated!" << std::endl; return false; }; // 行列生成 if (!gen_mtx(h, w, mtx)) { std::cout << "[ERROR] Matrix could not be generated!" << std::endl; return false; }; // 配列 v 計算 if (!calc_v(mtx, v)) { std::cout << "[ERROR] V array could not be calculated!" << std::endl; return false; }; // 配列 b 計算 if (!calc_b(v, b)) { std::cout << "[ERROR] B array could not be calculated!" << std::endl; return false; }; // 配列 a 計算 if (!calc_a(v, a)) { std::cout << "[ERROR] B array could not be calculated!" << std::endl; return false; }; // 配列 d 計算 if (!calc_d(v, d)) { std::cout << "[ERROR] D array could not be calculated!" << std::endl; return false; }; // 配列 c 計算 if (!calc_c(v, c)) { std::cout << "[ERROR] C array could not be calculated!" << std::endl; return false; }; // 補間 s = 0.1; for (i = 0; i <= (data_src[n - 1][0] - data_src[0][0]) / s; i++) { std::vector row; x = i * s; idx = search_i(x); y = a[idx] * std::pow(x - data_src[idx][0], 3.0) + b[idx] * std::pow(x - data_src[idx][0], 2.0) + c[idx] * (x - data_src[idx][0]) + d[idx]; row.push_back(x); row.push_back(y); data_res.push_back(row); } } catch (...) { return false; // 計算失敗 } return true; // 計算成功 } bool Calc::calc_a( std::vector& v, std::vector& a ) { std::size_t i; try { for (i = 0; i < n - 1; i++) a.push_back( (v[i + 1] - v[i]) / (6 * (data_src[i + 1][0] - data_src[i][0])) ); } catch (...) { return false; } return true; } bool Calc::calc_b( std::vector& v, std::vector& b ) { std::size_t i; try { for (i = 0; i < n - 1; i++) b.push_back(v[i] / 2.0); } catch (...) { return false; } return true; } bool Calc::calc_c( std::vector& v, std::vector& c ) { std::size_t i; try { for (i = 0; i < n - 1; i++) c.push_back( (data_src[i + 1][1] - data_src[i][1]) / (data_src[i + 1][0] - data_src[i][0]) - (data_src[i + 1][0] - data_src[i][0]) * (2.0 * v[i] + v[i + 1]) / 6.0 ); } catch (...) { return false; } return true; } bool Calc::calc_d( std::vector& v, std::vector& d ) { std::size_t i; try { for (i = 0; i < n; i++) d.push_back(data_src[i][1]); } catch (...) { return false; } return true; } bool Calc::calc_h( std::vector>& data, std::vector& h ) { std::size_t i; try { for (i = 0; i < n - 1; i++) h.push_back(data[i + 1][0] - data[i][0]); } catch (...) { return false; } return true; } bool Calc::calc_v( std::vector>& mtx, std::vector& v ) { std::size_t i; std::vector g; // Gauss-Jordan 計算用配列 try { if (!solve_gauss_jordan(mtx, g)) { std::cout << "[ERROR] Could not be solved by Gauss-Jordan!" << std::endl; return false; }; v.push_back(0.0); for (i = 0; i < n - 2; i++) v.push_back(g[i]); v.push_back(0.0); } catch (...) { return false; } return true; } bool Calc::calc_w( std::vector>& data, std::vector& h, std::vector& w ) { std::size_t i; try { for (i = 1; i < n - 1; i++) w.push_back( 6.0 * ((data[i + 1][1] - data[i][1]) / h[i] - (data[i][1] - data[i - 1][1]) / h[i - 1]) ) ; } catch (...) { return false; } return true; } bool Calc::gen_mtx( std::vector& h, std::vector& w, std::vector>& mtx ) { std::size_t i; std::size_t j; try { for (i = 0; i < n - 2; i++) { std::vector row; for (j = 0; j < n - 1; j++) row.push_back(0.0); mtx.push_back(row); } for (i = 0; i < n - 2; i++) { mtx[i][i] = 2.0 * (h[i] + h[i + 1]); mtx[i][n - 2] = w[i]; if (i == 0) continue; mtx[i - 1][i] = h[i]; mtx[i][i - 1] = h[i]; } } catch (...) { return false; } return true; } bool Calc::solve_gauss_jordan( std::vector>& mtx, std::vector& g ) { std::size_t i; std::size_t j; std::size_t k; double p; double d; try { for (k = 0; k < n - 2; k++) { p = mtx[k][k]; for (j = k; j < n - 1; j++) mtx[k][j] /= p; for (i = 0; i < n - 2; i++) { if (i == k) continue; d = mtx[i][k]; for (j = k; j < n - 1; j++) mtx[i][j] -= d * mtx[k][j]; } } for (i = 0; i < n - 2; i++) g.push_back(mtx[i][n - 2]); } catch (...) { return false; } return true; } unsigned int Calc::search_i(double x) { std::size_t i; std::size_t j; double k; try { i = 0.0; j = n - 1; while (i < j) { k = (i + j) / 2.0; if (data_src[k][0] < x) { i = k + 1; } else { j = k; } } if (i > 0.0) i--; return i; } catch (...) { return 0; } }
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 #ifndef SPLINE_INTERPOLATION_CALC_HPP_ #define SPLINE_INTERPOLATION_CALC_HPP_ #include class Calc { std::vector> data_src; // 元データ std::size_t n; // 与えられた点の数 bool calc_a( std::vector&, std::vector& ); // 配列 a 計算 bool calc_b( std::vector&, std::vector& ); // 配列 b 計算 bool calc_c( std::vector&, std::vector& ); // 配列 c 計算 bool calc_d( std::vector&, std::vector& ); // 配列 d 計算 bool calc_h( std::vector>&, std::vector& ); // 配列 h 計算 bool calc_v( std::vector>&, std::vector& ); // 配列 v 計算 bool calc_w( std::vector>&, std::vector&, std::vector& ); // 配列 w 計算 bool gen_mtx( std::vector&, std::vector&, std::vector>& ); // 計算用行列生成 bool solve_gauss_jordan( std::vector>&, std::vector& ); // Gauss-Jordan 解法 unsigned int search_i(double); // 計算対象インデックス検索 public: Calc(std::vector>& data_src) : data_src(data_src), n(get_size()) {} // コンストラクタ unsigned int get_size(); // 行列サイズ取得 bool interpolate_spline( std::vector>&, std::vector>& ); // 3次スプライン補間 }; #endif
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 #include "file.hpp" #include #include #include #include bool File::get_text(std::vector>& data) { try { // ファイル OPEN std::ifstream ifs(f_data); if (!ifs.is_open()) return false; // 読み込み失敗 // ファイル READ std::string buf; // 1行分バッファ while (getline(ifs, buf)) { std::vector rec; // 1行分ベクタ std::istringstream iss(buf); // 文字列ストリーム std::string field; // 1列分文字列 // 1行分文字列を1行分ベクタに追加 double s; while (iss >> s) rec.push_back(s); // １行分ベクタを data ベクタに追加 if (rec.size() != 0) data.push_back(rec); } } catch (...) { std::cerr << "EXCEPTION!" << std::endl; return false; } return true; // 読み込み成功 }
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 #ifndef SPLINE_INTERPOLATION_FILE_HPP_ #define SPLINE_INTERPOLATION_FILE_HPP_ #include #include #include class File { std::string f_data; public: File(std::string f_data) : f_data(f_data) {} bool get_text(std::vector>&); }; #endif
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 /*********************************************************** 3次スプライン補間 DATE AUTHOR VERSION 2020.07.31 mk-mode.com 1.00 新規作成 Copyright(C) 2020 mk-mode.com All Rights Reserved. ***********************************************************/ #include "calc.hpp" #include "file.hpp" #include // for EXIT_XXXX #include // for setprecision #include #include #include int main(int argc, char* argv[]) { std::string f_data; // データファイル名 std::vector> data_src; // データ配列（与えられた点） std::vector> data_res; // データ配列（計算結果） std::size_t i; // loop インデックス try { // コマンドライン引数のチェック if (argc != 2) { std::cerr << "[ERROR] Number of arguments is wrong!\n" << "[USAGE] ./spline_interpolation " << std::endl; return EXIT_FAILURE; } // ファイル名取得 f_data = argv[1]; // データ取得 File file(f_data); if (!file.get_text(data_src)) { std::cout << "[ERROR] Failed to read the file!" << std::endl; return EXIT_FAILURE; } // 計算用オプジェクトのインスタンス化 Calc calc(data_src); // 計算（補間） if (!calc.interpolate_spline(data_src, data_res)) { std::cout << "[ERROR] Failed to decompose!" << std::endl; return EXIT_FAILURE; } // 結果出力 std::cout << std::setw(10) << std::right << "x" << std::setw(10) << std::right << "y" << std::endl; std::cout << std::fixed << std::setprecision(4); for (i = 0; i < data_res.size(); i++) { std::cout << std::setw(10) << std::right << data_res[i][0] << std::setw(10) << std::right << data_res[i][1] << std::endl; } } catch (...) { std::cerr << "EXCEPTION!" << std::endl; return EXIT_FAILURE; } return EXIT_SUCCESS; }

### MariaMarinari commented Nov 10, 2022

Hi komasaru,

I am looking for a 3D cubic spline code that I could use for a model I am running. I was wondering if this code fitted a line in 3D or a surface?