C++ source code to interpolate by 3D-Spline.

calc.cpp

#include "calc.hpp"

#include <cmath>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <vector>

/**
 * @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<std::vector<double>>& data_src,
  std::vector<std::vector<double>>& data_res
) {
  std::vector<double> a;                 // 配列 a
  std::vector<double> b;                 // 配列 b
  std::vector<double> c;                 // 配列 c
  std::vector<double> d;                 // 配列 a
  std::vector<double> h;                 // 配列 h
  std::vector<double> v;                 // 配列 v
  std::vector<double> w;                 // 配列 w
  std::vector<std::vector<double>> 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<double> 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<double>& v,
  std::vector<double>& 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<double>& v,
  std::vector<double>& 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<double>& v,
  std::vector<double>& 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<double>& v,
  std::vector<double>& 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<std::vector<double>>& data,
  std::vector<double>& 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<std::vector<double>>& mtx,
  std::vector<double>& v
) {
  std::size_t i;
  std::vector<double> 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<std::vector<double>>& data,
  std::vector<double>& h,
  std::vector<double>& 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<double>& h,
  std::vector<double>& w,
  std::vector<std::vector<double>>& mtx
) {
  std::size_t i;
  std::size_t j;

  try {
    for (i = 0; i < n - 2; i++) {
      std::vector<double> 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<std::vector<double>>& mtx,
  std::vector<double>& 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;
  }
}

calc.hpp

#ifndef SPLINE_INTERPOLATION_CALC_HPP_
#define SPLINE_INTERPOLATION_CALC_HPP_

#include <vector>

class Calc {
  std::vector<std::vector<double>> data_src;  // 元データ
  std::size_t n;                              // 与えられた点の数
  bool calc_a(
    std::vector<double>&,
    std::vector<double>&
  );                                          // 配列 a 計算
  bool calc_b(
    std::vector<double>&,
    std::vector<double>&
  );                                          // 配列 b 計算
  bool calc_c(
    std::vector<double>&,
    std::vector<double>&
  );                                          // 配列 c 計算
  bool calc_d(
    std::vector<double>&,
    std::vector<double>&
  );                                          // 配列 d 計算
  bool calc_h(
    std::vector<std::vector<double>>&,
    std::vector<double>&
  );                                          // 配列 h 計算
  bool calc_v(
    std::vector<std::vector<double>>&,
    std::vector<double>&
  );                                          // 配列 v 計算
  bool calc_w(
    std::vector<std::vector<double>>&,
    std::vector<double>&,
    std::vector<double>&
  );                                          // 配列 w 計算
  bool gen_mtx(
    std::vector<double>&,
    std::vector<double>&,
    std::vector<std::vector<double>>&
  );                                          // 計算用行列生成
  bool solve_gauss_jordan(
    std::vector<std::vector<double>>&,
    std::vector<double>&
  );                                          // Gauss-Jordan 解法
  unsigned int search_i(double);              // 計算対象インデックス検索

public:
  Calc(std::vector<std::vector<double>>& data_src)
    : data_src(data_src), n(get_size()) {}    // コンストラクタ
  unsigned int get_size();                    // 行列サイズ取得
  bool interpolate_spline(
    std::vector<std::vector<double>>&,
    std::vector<std::vector<double>>&
  );                                          // 3次スプライン補間
};

#endif

file.cpp

#include "file.hpp"

#include <iostream>
#include <sstream>
#include <string>
#include <vector>

bool File::get_text(std::vector<std::vector<double>>& 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<double> 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;  // 読み込み成功
}

file.hpp

#ifndef SPLINE_INTERPOLATION_FILE_HPP_
#define SPLINE_INTERPOLATION_FILE_HPP_

#include <fstream>
#include <string>
#include <vector>

class File {
  std::string f_data;

public:
  File(std::string f_data) : f_data(f_data) {}
  bool get_text(std::vector<std::vector<double>>&);
};

#endif

spline_interpolation.cpp

/***********************************************************
  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 <cstdlib>   // for EXIT_XXXX
#include <iomanip>   // for setprecision
#include <iostream>
#include <string>
#include <vector>

int main(int argc, char* argv[]) {
  std::string f_data;                         // データファイル名
  std::vector<std::vector<double>> data_src;  // データ配列（与えられた点）
  std::vector<std::vector<double>> data_res;  // データ配列（計算結果）
  std::size_t i;                              // loop インデックス

  try {
    // コマンドライン引数のチェック
    if (argc != 2) {
      std::cerr << "[ERROR] Number of arguments is wrong!\n"
                << "[USAGE] ./spline_interpolation <file_name>"
                << 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;
}

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?