sasekazu

// jacobi法による実対称行列の対角化
// numeric.jsが必要
// 他の関数との依存性排除ver
// size: 行列サイズ
// intputA: 対角化される入力行列  ※対称行列でなければならない
// 返り値 {DiagonalMatrix:A, OrthodonalMatrix:P};
// A(入力時) = P A(対角化後) P.transpose が成り立つ
function jacobiDiagonalization(size, inputA, maxIteration, threshold) {

	var A=numeric.clone(inputA);

sasekazu

// numeric.jsを使用
// データ点を配列で渡す 例: points = [[x0,y0],[x1,y1],[x2,y2],[x3,y3]];
// y=c[0]*x+c[1] で近似し係数cを配列で返す
function leastSquare(points) {
	var A=numeric.rep([2, 2], 0);
	var b=numeric.rep([2], 0);
	for(var i=0; i<points.length; ++i) {
		A[0][0]+=points[i][0]*points[i][0];
		A[0][1]+=points[i][0];
		A[1][0]+=points[i][0];

sasekazu

// numeric.jsを使用
// データ点を配列で渡す 例: points = [[x0,y0],[x1,y1],[x2,y2],[x3,y3]];
// 多項式の係数を配列で返す
// dim=2 -> 直線
// dim=3 -> 二次曲線
function leastSquareN(points, dim) {
	if(points.length<dim) {
		return numeric.rep([dim], 0);
	}
	var A=numeric.rep([dim, dim], 0);

sasekazu

// ラグランジュ補間（と等価な多項式補間）
// 入力 N個の点列 points = [[x0,y0],[x1,x2],...,[xn-1,yn-1]]
// 出力 多項式の係数配列c
//     y = c[0]*x^N-1 + c[1]*x^N-2 + ... + c[N-2]*x + c[N-1]
function lagrangeInterpolation(points) {
	var dim=points.length;
	var A=numeric.rep([dim, dim], 0);
	var b=numeric.rep([dim], 0);
	for(var i=0; i<dim; ++i) {
		for(var j=0; j<dim; ++j) {

sasekazu

// 3次スプライン補間
// 参考：高橋大輔，数値計算，岩波書店，1996.
// 入力 N+1個の点列 points = [[x0,y0],[x1,x2],...,[xn,yn]]
// 出力 各データ点間の区間における多項式の係数を含む2次元配列
// x：x座標値で昇順ソートされた点列のx座標
// y：x座標値で昇順ソートされた点列のy座標
// a：区間[x_i+1 - x_i]における3次式の係数 配列サイズN-1
// b：区間[x_i+1 - x_i]における3次式の係数 配列サイズN-1
// c：区間[x_i+1 - x_i]における3次式の係数 配列サイズN-1
// d：区間[x_i+1 - x_i]における3次式の係数 配列サイズN-1

sasekazu

// 入力誤差モデルを用いた最尤推定による直線当てはめ
// 入力 データ点 points = [[x0,y0],[x1,y1],...,[xn-1,yn-1]]
// y = ax + b の a, b を格納した配列 [a,b] を返す
function lineFittingInputErrorModel(points) {

	// 点の数
	var N=points.length;	

	// xの重心

sasekazu

// サイズ NxN の行列 A の逆行列を計算する
// ピボット選択なし
// 実際にはnumeric javascript の　numeric.inverse() や numeric.solve() を使うべき
function inverseGaussJordan(A, N) {
	var U=numeric.clone(A);
	var Ainv=identityMatrix(N);
	var alpha;
	var diagValue;
	for(var i=0; i<N; ++i) {
		diagValue=U[i][i];

sasekazu

// ドロネー三角形分割関数
// ※簡潔だが効率の悪い実装
// 引数1 inputPoints: 入力点の座標 [[x1,y1],[x2,y2],....]
// 引数2,3,4,5 ymax, ymin, xmax, xmin: 入力点が含まれる領域の最大・最小座標
// 返り値：triangleのコネクティビティ 
// 返り値の例: [[0,1,2], [0,1,3], [1,2,4]..] の場合，
// inputPoints[0], inputPoints[1], inputPoints[2] が三角形1を構成する。…以下、同様
function DelaunayTriangulation(inputPoints, ymax, ymin, xmax, xmin) {

	var pos = numeric.clone(inputPoints);	// 点の数 x 2(x,y)

sasekazu

function fileSave(text, filename){
	var blob = new Blob([text],{"type" : "text/html"});
	var newElement = document.createElement("a");   // 新しいDOMを作成
	newElement.textContent= "download";             // テキストはdownload
	newElement.setAttribute('href', window.URL.createObjectURL(blob));  
	newElement.setAttribute('download', filename);
	document.body.appendChild(newElement);          // body直下に作成
}

sasekazu

#include <iostream>
#include <GL/glew.h>
#include <GL/glut.h>
#include <oglplus/all.hpp>

oglplus::Context* gl;

void Display(void);

int main(int argc, char* argv[]){