+/* \r
+ * PROJECT: NyARToolkit\r
+ * --------------------------------------------------------------------------------\r
+ * This work is based on the original ARToolKit developed by\r
+ * Hirokazu Kato\r
+ * Mark Billinghurst\r
+ * HITLab, University of Washington, Seattle\r
+ * http://www.hitl.washington.edu/artoolkit/\r
+ *\r
+ * The NyARToolkit is Java version ARToolkit class library.\r
+ * Copyright (C)2008 R.Iizuka\r
+ *\r
+ * This program is free software; you can redistribute it and/or\r
+ * modify it under the terms of the GNU General Public License\r
+ * as published by the Free Software Foundation; either version 2\r
+ * of the License, or (at your option) any later version.\r
+ * \r
+ * This program is distributed in the hope that it will be useful,\r
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of\r
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\r
+ * GNU General Public License for more details.\r
+ * \r
+ * You should have received a copy of the GNU General Public License\r
+ * along with this framework; if not, write to the Free Software\r
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
+ * \r
+ * For further information please contact.\r
+ * http://nyatla.jp/nyatoolkit/\r
+ * <airmail(at)ebony.plala.or.jp>\r
+ * \r
+ */\r
+package jp.nyatla.nyartoolkit.core;\r
+\r
+import jp.nyatla.nyartoolkit.NyARException;\r
+\r
+\r
+\r
+/**\r
+ * ARMat構造体に対応するクラス typedef struct { double *m; int row; int clm; }ARMat;\r
+ * \r
+ */\r
+public class NyARMat\r
+{\r
+ /**\r
+ * 配列サイズと行列サイズは必ずしも一致しないことに注意 返された配列のサイズを行列の大きさとして使わないこと!\r
+ * \r
+ */\r
+ protected double[][] m;\r
+\r
+ protected int clm, row;\r
+\r
+ /**\r
+ * デフォルトコンストラクタは機能しません。\r
+ * \r
+ * @throws NyARException\r
+ */\r
+ protected NyARMat() throws NyARException\r
+ {\r
+ throw new NyARException();\r
+ }\r
+\r
+ public NyARMat(int i_row, int i_clm)\r
+ {\r
+ m = new double[i_row][i_clm];\r
+ clm = i_clm;\r
+ row = i_row;\r
+ }\r
+\r
+ /**\r
+ * i_row x i_clmサイズの行列を格納できるように行列サイズを変更します。 実行後、行列の各値は不定になります。\r
+ * \r
+ * @param i_row\r
+ * @param i_clm\r
+ */\r
+ public void realloc(int i_row, int i_clm)\r
+ {\r
+ if (i_row <= this.m.length && i_clm <= this.m[0].length) {\r
+ // 十分な配列があれば何もしない。\r
+ } else {\r
+ // 不十分なら取り直す。\r
+ m = new double[i_row][i_clm];\r
+ }\r
+ this.clm = i_clm;\r
+ this.row = i_row;\r
+ }\r
+\r
+ public int getClm()\r
+ {\r
+ return clm;\r
+ }\r
+\r
+ public int getRow()\r
+ {\r
+ return row;\r
+ }\r
+\r
+ /**\r
+ * 行列をゼロクリアする。\r
+ */\r
+ public void zeroClear()\r
+ {\r
+ int i, i2;\r
+ // For順変更OK\r
+ for (i = row - 1; i >= 0; i--) {\r
+ for (i2 = clm - 1; i2 >= 0; i2--) {\r
+ m[i][i2] = 0.0;\r
+ }\r
+ }\r
+ }\r
+\r
+ /**\r
+ * i_copy_fromの内容を自分自身にコピーします。 高さ・幅は同一で無いと失敗します。\r
+ * \r
+ * @param i_copy_from\r
+ */\r
+ public void copyFrom(NyARMat i_copy_from) throws NyARException\r
+ {\r
+ // サイズ確認\r
+ if (this.row != i_copy_from.row || this.clm != i_copy_from.clm) {\r
+ throw new NyARException();\r
+ }\r
+ // 値コピー\r
+ for (int r = this.row - 1; r >= 0; r--) {\r
+ for (int c = this.clm - 1; c >= 0; c--) {\r
+ this.m[r][c] = i_copy_from.m[r][c];\r
+ }\r
+ }\r
+ }\r
+\r
+ public double[][] getArray()\r
+ {\r
+ return m;\r
+ }\r
+\r
+ // public void getRowVec(int i_row,NyARVec o_vec)\r
+ // {\r
+ // o_vec.set(this.m[i_row],this.clm);\r
+ // }\r
+ /**\r
+ * aとbの積を自分自身に格納する。arMatrixMul()の代替品\r
+ * \r
+ * @param a\r
+ * @param b\r
+ * @throws NyARException\r
+ */\r
+ public void matrixMul(NyARMat a, NyARMat b) throws NyARException\r
+ {\r
+ if (a.clm != b.row || this.row != a.row || this.clm != b.clm) {\r
+ throw new NyARException();\r
+ }\r
+ double w;\r
+ int r, c, i;\r
+ double[][] am = a.m, bm = b.m, dm = this.m;\r
+ // For順変更禁止\r
+ for (r = 0; r < this.row; r++) {\r
+ for (c = 0; c < this.clm; c++) {\r
+ w = 0.0;// dest.setARELEM0(r, c,0.0);\r
+ for (i = 0; i < a.clm; i++) {\r
+ w += am[r][i] * bm[i][c];// ARELEM0(dest, r, c) +=ARELEM0(a, r, i) * ARELEM0(b,i, c);\r
+ }\r
+ dm[r][c] = w;\r
+ }\r
+ }\r
+ }\r
+\r
+ private int[] wk_nos_matrixSelfInv = new int[50];\r
+\r
+ // private final static double matrixSelfInv_epsl=1.0e-10;\r
+ /**\r
+ * i_targetを逆行列に変換する。arMatrixSelfInv()と、arMatrixSelfInv_minv()関数を合成してあります。\r
+ * OPTIMIZE STEP[485->422]\r
+ * \r
+ * @param i_target\r
+ * 逆行列にする行列\r
+ * @return 逆行列があればTRUE/無ければFALSE\r
+ * \r
+ * @throws NyARException\r
+ */\r
+ public boolean matrixSelfInv() throws NyARException\r
+ {\r
+ double[][] ap = this.m;\r
+ int dimen = this.row;\r
+ int dimen_1 = dimen - 1;\r
+ double[] ap_n, ap_ip, ap_i;// wap;\r
+ int j, ip, nwork;\r
+ int[] nos = wk_nos_matrixSelfInv;// この関数で初期化される。\r
+ // double epsl;\r
+ double p, pbuf, work;\r
+\r
+ /* check size */\r
+ switch (dimen) {\r
+ case 0:\r
+ throw new NyARException();\r
+ case 1:\r
+ ap[0][0] = 1.0 / ap[0][0];// *ap = 1.0 / (*ap);\r
+ return true;/* 1 dimension */\r
+ }\r
+\r
+ for (int n = 0; n < dimen; n++) {\r
+ nos[n] = n;\r
+ }\r
+\r
+ /*\r
+ * nyatla memo ipが定まらないで計算が行われる場合があるので挿入。 ループ内で0初期化していいかが判らない。\r
+ */\r
+ ip = 0;\r
+ // For順変更禁止\r
+ for (int n = 0; n < dimen; n++) {\r
+ ap_n = ap[n];// wcp = ap + n * rowa;\r
+ p = 0.0;\r
+ for (int i = n; i < dimen; i++) {// for(i = n, wap = wcp, p =\r
+ // 0.0; i < dimen ; i++, wap +=\r
+ // rowa)\r
+ if (p < (pbuf = Math.abs(ap[i][0]))) {\r
+ p = pbuf;\r
+ ip = i;\r
+ }\r
+ }\r
+ // if (p <= matrixSelfInv_epsl){\r
+ if (p == 0.0) {\r
+ return false;\r
+ // throw new NyARException();\r
+ }\r
+\r
+ nwork = nos[ip];\r
+ nos[ip] = nos[n];\r
+ nos[n] = nwork;\r
+\r
+ ap_ip = ap[ip];\r
+ for (j = 0; j < dimen; j++) {// for(j = 0, wap = ap + ip * rowa,\r
+ // wbp = wcp; j < dimen ; j++) {\r
+ work = ap_ip[j]; // work = *wap;\r
+ ap_ip[j] = ap_n[j];\r
+ ap_n[j] = work;\r
+ }\r
+\r
+ work = ap_n[0];\r
+ for (j = 0; j < dimen_1; j++) {// for(j = 1, wap = wcp, work =\r
+ // *wcp; j < dimen ; j++, wap++)\r
+ ap_n[j] = ap_n[j + 1] / work;// *wap = *(wap + 1) / work;\r
+ }\r
+ ap_n[j] = 1.0 / work;// *wap = 1.0 / work;\r
+ for (int i = 0; i < dimen; i++) {\r
+ if (i != n) {\r
+ ap_i = ap[i];// wap = ap + i * rowa;\r
+\r
+ work = ap_i[0];\r
+ for (j = 0; j < dimen_1; j++) {// for(j = 1, wbp = wcp,work = *wap;j < dimen ;j++, wap++, wbp++)\r
+ ap_i[j] = ap_i[j + 1] - work * ap_n[j];// wap = *(wap +1) - work *(*wbp);\r
+ }\r
+ ap_i[j] = -work * ap_n[j];// *wap = -work * (*wbp);\r
+ }\r
+ }\r
+ }\r
+\r
+ for (int n = 0; n < dimen; n++) {\r
+ for (j = n; j < dimen; j++) {\r
+ if (nos[j] == n) {\r
+ break;\r
+ }\r
+ }\r
+ nos[j] = nos[n];\r
+ for (int i = 0; i < dimen; i++) {// for(i = 0, wap = ap + j, wbp\r
+ // = ap + n; i < dimen ;i++, wap\r
+ // += rowa, wbp += rowa) {\r
+ ap_i = ap[i];\r
+ work = ap_i[j];// work = *wap;\r
+ ap_i[j] = ap_i[n];// *wap = *wbp;\r
+ ap_i[n] = work;// *wbp = work;\r
+ }\r
+ }\r
+ return true;\r
+ }\r
+\r
+ /**\r
+ * sourceの転置行列をdestに得る。arMatrixTrans()の代替品\r
+ * \r
+ * @param dest\r
+ * @param source\r
+ * @return\r
+ */\r
+ public static void matrixTrans(NyARMat dest, NyARMat source) throws NyARException\r
+ {\r
+ if (dest.row != source.clm || dest.clm != source.row) {\r
+ throw new NyARException();\r
+ }\r
+ NyARException.trap("未チェックのパス");\r
+ // For順変更禁止\r
+ for (int r = 0; r < dest.row; r++) {\r
+ for (int c = 0; c < dest.clm; c++) {\r
+ dest.m[r][c] = source.m[c][r];\r
+ }\r
+ }\r
+ }\r
+\r
+ /**\r
+ * unitを単位行列に初期化する。arMatrixUnitの代替品\r
+ * \r
+ * @param unit\r
+ */\r
+ public static void matrixUnit(NyARMat unit) throws NyARException\r
+ {\r
+ if (unit.row != unit.clm) {\r
+ throw new NyARException();\r
+ }\r
+ NyARException.trap("未チェックのパス");\r
+ // For順変更禁止\r
+ for (int r = 0; r < unit.getRow(); r++) {\r
+ for (int c = 0; c < unit.getClm(); c++) {\r
+ if (r == c) {\r
+ unit.m[r][c] = 1.0;\r
+ } else {\r
+ unit.m[r][c] = 0.0;\r
+ }\r
+ }\r
+ }\r
+ }\r
+\r
+ /**\r
+ * sourceの内容を自身に複製する。 Optimized 2008.04.19\r
+ * \r
+ * @param i_source\r
+ * @return\r
+ */\r
+ public void matrixDup(NyARMat i_source) throws NyARException\r
+ {\r
+ // 自身の配列サイズを相手のそれより大きいことを保障する。\r
+ this.realloc(i_source.row, i_source.clm);\r
+ // 内容を転写\r
+ int r, c;\r
+ double[][] src_m, dest_m;\r
+ src_m = i_source.m;\r
+ dest_m = this.m;\r
+ // コピーはFor順を変えてもOK\r
+ for (r = this.row - 1; r >= 0; r--) {\r
+ for (c = this.clm - 1; c >= 0; c--) {\r
+ dest_m[r][c] = src_m[r][c];\r
+ }\r
+ }\r
+ }\r
+\r
+ public NyARMat matrixAllocDup() throws NyARException\r
+ {\r
+ NyARMat result = new NyARMat(this.row, this.clm);\r
+ // コピー\r
+ int r, c;\r
+ double[][] dest_m, src_m;\r
+ dest_m = result.m;\r
+ src_m = this.m;\r
+ // コピーはFor順を変えてもOK\r
+ for (r = this.row - 1; r >= 0; r--) {\r
+ for (c = this.clm - 1; c >= 0; c--) {\r
+ dest_m[r][c] = src_m[r][c];\r
+ }\r
+ }\r
+ return result;\r
+ }\r
+\r
+ /**\r
+ * arMatrixInv関数の代替品です。 destにsourceの逆行列を返します。\r
+ * \r
+ * @param dest\r
+ * @param source\r
+ * @throws NyARException\r
+ */\r
+ public static void matrixInv(NyARMat dest, NyARMat source)\r
+ throws NyARException\r
+ {\r
+ NyARException.trap("未チェックのパス");\r
+ dest.matrixDup(source);\r
+\r
+ NyARException.trap("未チェックのパス");\r
+ dest.matrixSelfInv();\r
+ }\r
+\r
+ public NyARMat matrixAllocInv() throws NyARException\r
+ {\r
+ NyARException.trap("未チェックのパス");\r
+ NyARMat result = matrixAllocDup();\r
+\r
+ NyARException.trap("未チェックのパス");\r
+ result.matrixSelfInv();\r
+ return result;\r
+ }\r
+\r
+ /**\r
+ * dim x dim の単位行列を作る。\r
+ * \r
+ * @param dim\r
+ * @return\r
+ * @throws NyARException\r
+ */\r
+ public static NyARMat matrixAllocUnit(int dim) throws NyARException\r
+ {\r
+ NyARException.trap("未チェックのパス");\r
+ NyARMat result = new NyARMat(dim, dim);\r
+ NyARException.trap("未チェックのパス");\r
+ NyARMat.matrixUnit(result);\r
+ return result;\r
+ }\r
+\r
+ /**\r
+ * arMatrixDispの代替品\r
+ * \r
+ * @param m\r
+ * @return\r
+ */\r
+ public int matrixDisp() throws NyARException\r
+ {\r
+ NyARException.trap("未チェックのパス");\r
+ System.out.println(" === matrix (" + row + "," + clm + ") ===");// printf(" ===matrix (%d,%d) ===\n", m->row, m->clm);\r
+ for (int r = 0; r < row; r++) {// for(int r = 0; r < m->row; r++) {\r
+ System.out.print(" |");// printf(" |");\r
+ for (int c = 0; c < clm; c++) {// for(int c = 0; c < m->clm; c++) {\r
+ System.out.print(" " + m[r][c]);// printf(" %10g", ARELEM0(m, r, c));\r
+ }\r
+ System.out.println(" |");// printf(" |\n");\r
+ }\r
+ System.out.println(" ======================");// printf(" ======================\n");\r
+ return 0;\r
+ }\r
+\r
+ private static final double PCA_EPS = 1e-6; // #define EPS 1e-6\r
+\r
+ private static final int PCA_MAX_ITER = 100; // #define MAX_ITER 100\r
+\r
+ private static final double PCA_VZERO = 1e-16; // #define VZERO 1e-16\r
+\r
+ /**\r
+ * static int EX( ARMat *input, ARVec *mean )の代替関数 Optimize:STEP:[144->110]\r
+ * \r
+ * @param input\r
+ * @param mean\r
+ * @return\r
+ * @throws NyARException\r
+ */\r
+ private void PCA_EX(NyARVec mean) throws NyARException\r
+ {\r
+ int lrow, lclm;\r
+ int i, i2;\r
+ lrow = this.row;\r
+ lclm = this.clm;\r
+ double[][] lm = this.m;\r
+\r
+ if (lrow <= 0 || lclm <= 0) {\r
+ throw new NyARException();\r
+ }\r
+ if (mean.getClm() != lclm) {\r
+ throw new NyARException();\r
+ }\r
+ // double[] mean_array=mean.getArray();\r
+ // mean.zeroClear();\r
+ final double[] mean_array = mean.getArray();\r
+ double w;\r
+ // For順変更禁止\r
+ for (i2 = 0; i2 < lclm; i2++) {\r
+ w = 0.0;\r
+ for (i = 0; i < lrow; i++) {\r
+ // *(v++) += *(m++);\r
+ w += lm[i][i2];\r
+ }\r
+ mean_array[i2] = w / lrow;// mean->v[i] /= row;\r
+ }\r
+ }\r
+\r
+ /**\r
+ * static int CENTER( ARMat *inout, ARVec *mean )の代替関数\r
+ * \r
+ * @param inout\r
+ * @param mean\r
+ * @return\r
+ */\r
+ private static void PCA_CENTER(NyARMat inout, NyARVec mean)throws NyARException\r
+ {\r
+ double[] v;\r
+ int row, clm;\r
+\r
+ row = inout.getRow();\r
+ clm = inout.getClm();\r
+ if (mean.getClm() != clm) {\r
+ throw new NyARException();\r
+ }\r
+ double[][] im = inout.m;\r
+ double[] im_i;\r
+ double w0, w1;\r
+ v = mean.getArray();\r
+ // 特にパフォーマンスが劣化するclm=1と2ときだけ、別パスで処理します。\r
+ switch (clm) {\r
+ case 1:\r
+ w0 = v[0];\r
+ for (int i = 0; i < row; i++) {\r
+ im[i][0] -= w0;\r
+ }\r
+ break;\r
+ case 2:\r
+ w0 = v[0];\r
+ w1 = v[1];\r
+ for (int i = 0; i < row; i++) {\r
+ im_i = im[i];\r
+ im_i[0] -= w0;\r
+ im_i[1] -= w1;\r
+ }\r
+ break;\r
+ default:\r
+ for (int i = 0; i < row; i++) {\r
+ im_i = im[i];\r
+ for (int j = 0; j < clm; j++) {\r
+ // *(m++) -= *(v++);\r
+ im_i[j] -= v[j];\r
+ }\r
+ }\r
+ break;\r
+ }\r
+ return;\r
+ }\r
+\r
+ /**\r
+ * int x_by_xt( ARMat *input, ARMat *output )の代替関数\r
+ * \r
+ * @param input\r
+ * @param output\r
+ * @throws NyARException\r
+ */\r
+ private static void PCA_x_by_xt(NyARMat input, NyARMat output) throws NyARException\r
+ {\r
+ NyARException.trap("動作未チェック/配列化未チェック");\r
+ int row, clm;\r
+ // double[][] out;\r
+ double[] in1, in2;\r
+\r
+ NyARException.trap("未チェックのパス");\r
+ row = input.row;\r
+ clm = input.clm;\r
+ NyARException.trap("未チェックのパス");\r
+ if (output.row != row || output.clm != row) {\r
+ throw new NyARException();\r
+ }\r
+\r
+ // out = output.getArray();\r
+ for (int i = 0; i < row; i++) {\r
+ for (int j = 0; j < row; j++) {\r
+ if (j < i) {\r
+ NyARException.trap("未チェックのパス");\r
+ output.m[i][j] = output.m[j][i];// *out =\r
+ // output->m[j*row+i];\r
+ } else {\r
+ NyARException.trap("未チェックのパス");\r
+ in1 = input.m[i];// input.getRowArray(i);//in1 = &(input->m[clm*i]);\r
+ in2 = input.m[j];// input.getRowArray(j);//in2 = &(input->m[clm*j]);\r
+ output.m[i][j] = 0;// *out = 0.0;\r
+ for (int k = 0; k < clm; k++) {\r
+ output.m[i][j] += (in1[k] * in2[k]);// *out += *(in1++)\r
+ // * *(in2++);\r
+ }\r
+ }\r
+ // out.incPtr();\r
+ }\r
+ }\r
+ }\r
+\r
+ /**\r
+ * static int xt_by_x( ARMat *input, ARMat *output )の代替関数\r
+ * Optimize:2008.04.19\r
+ * \r
+ * @param input\r
+ * @param i_output\r
+ * @throws NyARException\r
+ */\r
+ private static void PCA_xt_by_x(NyARMat input, NyARMat i_output) throws NyARException\r
+ {\r
+ double[] in_;\r
+ int row, clm;\r
+\r
+ row = input.row;\r
+ clm = input.clm;\r
+ if (i_output.row != clm || i_output.clm != clm) {\r
+ throw new NyARException();\r
+ }\r
+\r
+ int k, j;\r
+ double[][] out_m = i_output.m;\r
+ double w;\r
+ for (int i = 0; i < clm; i++) {\r
+ for (j = 0; j < clm; j++) {\r
+ if (j < i) {\r
+ out_m[i][j] = out_m[j][i];// *out = output->m[j*clm+i];\r
+ } else {\r
+ w = 0.0;// *out = 0.0;\r
+ for (k = 0; k < row; k++) {\r
+ in_ = input.m[k];// in=input.getRowArray(k);\r
+ w += (in_[i] * in_[j]);// *out += *in1 * *in2;\r
+ }\r
+ out_m[i][j] = w;\r
+ }\r
+ }\r
+ }\r
+ }\r
+\r
+ private final NyARVec wk_PCA_QRM_ev = new NyARVec(1);\r
+\r
+ /**\r
+ * static int QRM( ARMat *a, ARVec *dv )の代替関数\r
+ * \r
+ * @param a\r
+ * @param dv\r
+ * @throws NyARException\r
+ */\r
+ private void PCA_QRM(NyARVec dv) throws NyARException\r
+ {\r
+ double w, t, s, x, y, c;\r
+ int dim, iter;\r
+ double[] dv_array = dv.getArray();\r
+\r
+ dim = this.row;\r
+ if (dim != this.clm || dim < 2) {\r
+ throw new NyARException();\r
+ }\r
+ if (dv.getClm() != dim) {\r
+ throw new NyARException();\r
+ }\r
+\r
+ NyARVec ev = this.wk_PCA_QRM_ev;\r
+ ev.realloc(dim);\r
+ double[] ev_array = ev.getArray();\r
+ if (ev == null) {\r
+ throw new NyARException();\r
+ }\r
+ final double[][] L_m = this.m;\r
+ this.vecTridiagonalize(dv, ev, 1);\r
+\r
+ ev_array[0] = 0.0;// ev->v[0] = 0.0;\r
+ for (int h = dim - 1; h > 0; h--) {\r
+ int j = h;\r
+ while (j > 0&& Math.abs(ev_array[j]) > PCA_EPS* (Math.abs(dv_array[j - 1]) + Math.abs(dv_array[j]))) {// while(j>0 && fabs(ev->v[j]) >EPS*(fabs(dv->v[j-1])+fabs(dv->v[j])))\r
+ // j--;\r
+ j--;\r
+ }\r
+ if (j == h) {\r
+ continue;\r
+ }\r
+ iter = 0;\r
+ do {\r
+ iter++;\r
+ if (iter > PCA_MAX_ITER) {\r
+ break;\r
+ }\r
+ w = (dv_array[h - 1] - dv_array[h]) / 2;// w = (dv->v[h-1] -dv->v[h]) / 2;//ここ?\r
+ t = ev_array[h] * ev_array[h];// t = ev->v[h] * ev->v[h];\r
+ s = Math.sqrt(w * w + t);\r
+ if (w < 0) {\r
+ s = -s;\r
+ }\r
+ x = dv_array[j] - dv_array[h] + t / (w + s);// x = dv->v[j] -dv->v[h] +t/(w+s);\r
+ y = ev_array[j + 1];// y = ev->v[j+1];\r
+ for (int k = j; k < h; k++) {\r
+ if (Math.abs(x) >= Math.abs(y)) {\r
+ if (Math.abs(x) > PCA_VZERO) {\r
+ t = -y / x;\r
+ c = 1 / Math.sqrt(t * t + 1);\r
+ s = t * c;\r
+ } else {\r
+ c = 1.0;\r
+ s = 0.0;\r
+ }\r
+ } else {\r
+ t = -x / y;\r
+ s = 1.0 / Math.sqrt(t * t + 1);\r
+ c = t * s;\r
+ }\r
+ w = dv_array[k] - dv_array[k + 1];// w = dv->v[k] -dv->v[k+1];\r
+ t = (w * s + 2 * c * ev_array[k + 1]) * s;// t = (w * s +2 * c *ev->v[k+1]) *s;\r
+ dv_array[k] -= t;// dv->v[k] -= t;\r
+ dv_array[k + 1] += t;// dv->v[k+1] += t;\r
+ if (k > j) {\r
+ NyARException.trap("未チェックパス");\r
+ {\r
+ ev_array[k] = c * ev_array[k] - s * y;// ev->v[k]= c *ev->v[k]- s * y;\r
+ }\r
+ }\r
+ ev_array[k + 1] += s * (c * w - 2 * s * ev_array[k + 1]);// ev->v[k+1]+= s * (c* w- 2* s *ev->v[k+1]);\r
+\r
+ for (int i = 0; i < dim; i++) {\r
+ x = L_m[k][i];// x = a->m[k*dim+i];\r
+ y = L_m[k + 1][i];// y = a->m[(k+1)*dim+i];\r
+ L_m[k][i] = c * x - s * y;// a->m[k*dim+i] = c * x - s* y;\r
+ L_m[k + 1][i] = s * x + c * y;// a->m[(k+1)*dim+i] = s* x + c * y;\r
+ }\r
+ if (k < h - 1) {\r
+ NyARException.trap("未チェックパス");\r
+ {\r
+ x = ev_array[k + 1];// x = ev->v[k+1];\r
+ y = -s * ev_array[k + 2];// y = -s * ev->v[k+2];\r
+ ev_array[k + 2] *= c;// ev->v[k+2] *= c;\r
+ }\r
+ }\r
+ }\r
+ } while (Math.abs(ev_array[h]) > PCA_EPS\r
+ * (Math.abs(dv_array[h - 1]) + Math.abs(dv_array[h])));\r
+ }\r
+ for (int k = 0; k < dim - 1; k++) {\r
+ int h = k;\r
+ t = dv_array[h];// t = dv->v[h];\r
+ for (int i = k + 1; i < dim; i++) {\r
+ if (dv_array[i] > t) {// if( dv->v[i] > t ) {\r
+ h = i;\r
+ t = dv_array[h];// t = dv->v[h];\r
+ }\r
+ }\r
+ dv_array[h] = dv_array[k];// dv->v[h] = dv->v[k];\r
+ dv_array[k] = t;// dv->v[k] = t;\r
+ this.flipRow(h, k);\r
+ }\r
+ }\r
+\r
+ /**\r
+ * i_row_1番目の行と、i_row_2番目の行を入れ替える。\r
+ * \r
+ * @param i_row_1\r
+ * @param i_row_2\r
+ */\r
+ private void flipRow(int i_row_1, int i_row_2)\r
+ {\r
+ int i;\r
+ double w;\r
+ double[] r1 = this.m[i_row_1], r2 = this.m[i_row_2];\r
+ // For順変更OK\r
+ for (i = clm - 1; i >= 0; i--) {\r
+ w = r1[i];\r
+ r1[i] = r2[i];\r
+ r2[i] = w;\r
+ }\r
+ }\r
+\r
+ /**\r
+ * static int EV_create( ARMat *input, ARMat *u, ARMat *output, ARVec *ev\r
+ * )の代替関数\r
+ * \r
+ * @param input\r
+ * @param u\r
+ * @param output\r
+ * @param ev\r
+ * @throws NyARException\r
+ */\r
+ private static void PCA_EV_create(NyARMat input, NyARMat u, NyARMat output,NyARVec ev) throws NyARException\r
+ {\r
+ NyARException.trap("未チェックのパス");\r
+ int row, clm;\r
+ row = input.row;// row = input->row;\r
+ clm = input.clm;// clm = input->clm;\r
+ if (row <= 0 || clm <= 0) {\r
+ throw new NyARException();\r
+ }\r
+ if (u.row != row || u.clm != row) {// if( u->row != row || u->clm !=\r
+ // row ){\r
+ throw new NyARException();\r
+ }\r
+ if (output.row != row || output.clm != clm) {// if( output->row !=\r
+ // row || output->clm !=\r
+ // clm ){\r
+ throw new NyARException();\r
+ }\r
+ if (ev.getClm() != row) {// if( ev->clm != row ){\r
+ throw new NyARException();\r
+ }\r
+ double[][] m, in_;\r
+ double[] m1, ev_array;\r
+ double sum, work;\r
+\r
+ NyARException.trap("未チェックのパス");\r
+ m = output.m;// m = output->m;\r
+ in_ = input.m;\r
+ int i;\r
+ ev_array = ev.getArray();\r
+ for (i = 0; i < row; i++) {\r
+ NyARException.trap("未チェックのパス");\r
+ if (ev_array[i] < PCA_VZERO) {// if( ev->v[i] < VZERO ){\r
+ break;\r
+ }\r
+ NyARException.trap("未チェックのパス");\r
+ work = 1 / Math.sqrt(Math.abs(ev_array[i]));// work = 1 /\r
+ // sqrt(fabs(ev->v[i]));\r
+ for (int j = 0; j < clm; j++) {\r
+ sum = 0.0;\r
+ m1 = u.m[i];// m1 = &(u->m[i*row]);\r
+ // m2=input.getPointer(j);//m2 = &(input->m[j]);\r
+ for (int k = 0; k < row; k++) {\r
+ sum += m1[k] + in_[k][j];// sum += *m1 * *m2;\r
+ // m1.incPtr(); //m1++;\r
+ // m2.addPtr(clm);//m2 += clm;\r
+ }\r
+ m1[j] = sum * work;// *(m++) = sum * work;\r
+ // {//*(m++) = sum * work;\r
+ // m.set(sum * work);\r
+ // m.incPtr();}\r
+ }\r
+ }\r
+ for (; i < row; i++) {\r
+ NyARException.trap("未チェックのパス");\r
+ ev_array[i] = 0.0;// ev->v[i] = 0.0;\r
+ for (int j = 0; j < clm; j++) {\r
+ m[i][j] = 0.0;\r
+ // m.set(0.0);//*(m++) = 0.0;\r
+ // m.incPtr();\r
+ }\r
+ }\r
+ }\r
+\r
+ private NyARMat wk_PCA_PCA_u = null;\r
+\r
+ /**\r
+ * static int PCA( ARMat *input, ARMat *output, ARVec *ev )\r
+ * \r
+ * @param output\r
+ * @param o_ev\r
+ * @throws NyARException\r
+ */\r
+ private void PCA_PCA(NyARMat o_output, NyARVec o_ev) throws NyARException\r
+ {\r
+\r
+ int l_row, l_clm, min;\r
+ double[] ev_array = o_ev.getArray();\r
+\r
+ l_row = this.row;// row = input->row;\r
+ l_clm = this.clm;// clm = input->clm;\r
+ min = (l_clm < l_row) ? l_clm : l_row;\r
+ if (l_row < 2 || l_clm < 2) {\r
+ throw new NyARException();\r
+ }\r
+ if (o_output.clm != this.clm) {// if( output->clm != input->clm ){\r
+ throw new NyARException();\r
+ }\r
+ if (o_output.row != min) {// if( output->row != min ){\r
+ throw new NyARException();\r
+ }\r
+ if (o_ev.getClm() != min) {// if( ev->clm != min ){\r
+ throw new NyARException();\r
+ }\r
+\r
+ NyARMat u;// u =new NyARMat( min, min );\r
+ if (this.wk_PCA_PCA_u == null) {\r
+ u = new NyARMat(min, min);\r
+ this.wk_PCA_PCA_u = u;\r
+ } else {\r
+ u = this.wk_PCA_PCA_u;\r
+ u.realloc(min, min);\r
+ }\r
+\r
+ if (l_row < l_clm) {\r
+ NyARException.trap("未チェックのパス");\r
+ PCA_x_by_xt(this, u);// if(x_by_xt( input, u ) < 0 ) {\r
+ } else {\r
+ PCA_xt_by_x(this, u);// if(xt_by_x( input, u ) < 0 ) {\r
+ }\r
+ u.PCA_QRM(o_ev);\r
+\r
+ double[][] m1, m2;\r
+ if (l_row < l_clm) {\r
+ NyARException.trap("未チェックのパス");\r
+ PCA_EV_create(this, u, o_output, o_ev);\r
+ } else {\r
+ m1 = u.m;// m1 = u->m;\r
+ m2 = o_output.m;// m2 = output->m;\r
+ int i;\r
+ for (i = 0; i < min; i++) {\r
+ if (ev_array[i] < PCA_VZERO) {// if( ev->v[i] < VZERO ){\r
+ break;\r
+ }\r
+ for (int j = 0; j < min; j++) {\r
+ m2[i][j] = m1[i][j];// *(m2++) = *(m1++);\r
+ }\r
+ }\r
+ for (; i < min; i++) {// for( ; i < min; i++){\r
+ // コードを見た限りあってそうだからコメントアウト(2008/03/26)NyARException.trap("未チェックのパス");\r
+ ev_array[i] = 0.0;// ev->v[i] = 0.0;\r
+ for (int j = 0; j < min; j++) {\r
+ m2[i][j] = 0.0;// *(m2++) = 0.0;\r
+ }\r
+ }\r
+ }\r
+ }\r
+\r
+ private NyARMat wk_work_matrixPCA = null;\r
+\r
+ /**\r
+ * int arMatrixPCA( ARMat *input, ARMat *evec, ARVec *ev, ARVec *mean );\r
+ * 関数の置き換え。input引数がthisになる。 Optimize:2008.04.19\r
+ * \r
+ * @param o_evec\r
+ * @param o_ev\r
+ * \r
+ * @param mean\r
+ * @throws NyARException\r
+ */\r
+ public void matrixPCA(NyARMat o_evec, NyARVec o_ev, NyARVec mean)throws NyARException\r
+ {\r
+ double srow, sum;\r
+ int l_row, l_clm;\r
+ int check;\r
+\r
+ l_row = this.row;// row = input->row;\r
+ l_clm = this.clm;// clm = input->clm;\r
+ check = (l_row < l_clm) ? l_row : l_clm;\r
+ if (l_row < 2 || l_clm < 2) {\r
+ throw new NyARException();\r
+ }\r
+ if (o_evec.clm != l_clm || o_evec.row != check) {// if( evec->clm !=\r
+ // input->clm ||\r
+ // evec->row !=\r
+ // check ){\r
+ throw new NyARException();\r
+ }\r
+ if (o_ev.getClm() != check) {// if( ev->clm != check ){\r
+ throw new NyARException();\r
+ }\r
+ if (mean.getClm() != l_clm) {// if( mean->clm != input->clm ){\r
+ throw new NyARException();\r
+ }\r
+\r
+ // 自分の内容をワークにコピー(高速化の為に、1度作ったインスタンスは使いまわす)\r
+ NyARMat work;\r
+ if (this.wk_work_matrixPCA == null) {\r
+ work = this.matrixAllocDup();\r
+ this.wk_work_matrixPCA = work;\r
+ } else {\r
+ work = this.wk_work_matrixPCA;\r
+ work.matrixDup(this);// arMatrixAllocDup( input );work =\r
+ // arMatrixAllocDup( input );\r
+ }\r
+\r
+ srow = Math.sqrt((double) l_row);\r
+ work.PCA_EX(mean);\r
+\r
+ PCA_CENTER(work, mean);\r
+\r
+ int i, j;\r
+ // For順変更OK\r
+ for (i = 0; i < l_row; i++) {\r
+ for (j = 0; j < l_clm; j++) {\r
+ work.m[i][j] /= srow;// work->m[i] /= srow;\r
+ }\r
+ }\r
+\r
+ work.PCA_PCA(o_evec, o_ev);\r
+\r
+ sum = 0.0;\r
+ double[] ev_array = o_ev.getArray();\r
+ int ev_clm = o_ev.getClm();\r
+ // For順変更禁止\r
+ for (i = 0; i < ev_clm; i++) {// for(int i = 0; i < ev->clm; i++ ){\r
+ sum += ev_array[i];// sum += ev->v[i];\r
+ }\r
+ // For順変更禁止\r
+ for (i = 0; i < ev_clm; i++) {// for(int i = 0; i < ev->clm; i++ ){\r
+ ev_array[i] /= sum;// ev->v[i] /= sum;\r
+ }\r
+ }\r
+\r
+ /* int arMatrixPCA2( ARMat *input, ARMat *evec, ARVec *ev ); */\r
+ public static void arMatrixPCA2(NyARMat input, NyARMat evec, NyARVec ev) throws NyARException\r
+ {\r
+ NyARException.trap("未チェックのパス");\r
+ NyARMat work;\r
+ // double srow; // unreferenced\r
+ double sum;\r
+ int row, clm;\r
+ int check;\r
+\r
+ row = input.row;// row = input->row;\r
+ clm = input.clm;// clm = input->clm;\r
+ check = (row < clm) ? row : clm;\r
+ if (row < 2 || clm < 2) {\r
+ throw new NyARException();\r
+ }\r
+ if (evec.getClm() != input.clm || evec.row != check) {// if( evec->clm!= input->clm|| evec->row!= check ){\r
+ throw new NyARException();\r
+ }\r
+ if (ev.getClm() != check) {// if( ev->clm != check ){\r
+ throw new NyARException();\r
+ }\r
+\r
+ NyARException.trap("未チェックのパス");\r
+ work = input.matrixAllocDup();\r
+\r
+ NyARException.trap("未チェックパス");\r
+ work.PCA_PCA(evec, ev);// rval = PCA( work, evec, ev );\r
+ sum = 0.0;\r
+ double[] ev_array = ev.getArray();\r
+ for (int i = 0; i < ev.getClm(); i++) {// for( i = 0; i < ev->clm; i++\r
+ // ){\r
+ NyARException.trap("未チェックパス");\r
+ sum += ev_array[i];// sum += ev->v[i];\r
+ }\r
+ for (int i = 0; i < ev.getClm(); i++) {// for(int i = 0; i < ev->clm;i++ ){\r
+ NyARException.trap("未チェックパス");\r
+ ev_array[i] /= sum;// ev->v[i] /= sum;\r
+ }\r
+ return;\r
+ }\r
+\r
+ public static NyARMat matrixAllocMul(NyARMat a, NyARMat b) throws NyARException\r
+ {\r
+ NyARException.trap("未チェックのパス");\r
+ NyARMat dest = new NyARMat(a.row, b.clm);\r
+ NyARException.trap("未チェックのパス");\r
+ dest.matrixMul(a, b);\r
+ return dest;\r
+ }\r
+\r
+ /* static double mdet(double *ap, int dimen, int rowa) */\r
+ private static double Det_mdet(double[][] ap, int dimen, int rowa) throws NyARException\r
+ {\r
+ NyARException.trap("動作未チェック/配列化未チェック");\r
+ double det = 1.0;\r
+ double work;\r
+ int is_ = 0;\r
+ int mmax;\r
+\r
+ for (int k = 0; k < dimen - 1; k++) {\r
+ mmax = k;\r
+ for (int i = k + 1; i < dimen; i++) {\r
+ // if (Math.abs(arMatrixDet_MATRIX_get(ap, i, k, rowa)) >\r
+ // Math.abs(arMatrixDet_MATRIX_get(ap, mmax, k, rowa))){\r
+ if (Math.abs(ap[i][k]) > Math.abs(ap[mmax][k])) {\r
+ mmax = i;\r
+ }\r
+ }\r
+ if (mmax != k) {\r
+ for (int j = k; j < dimen; j++) {\r
+ work = ap[k][j];// work = MATRIX(ap, k, j, rowa);\r
+ ap[k][j] = ap[mmax][j];// MATRIX(ap, k, j, rowa) =MATRIX(ap, mmax, j, rowa);\r
+ ap[mmax][j] = work;// MATRIX(ap, mmax, j, rowa) = work;\r
+ }\r
+ is_++;\r
+ }\r
+ for (int i = k + 1; i < dimen; i++) {\r
+ work = ap[i][k] / ap[k][k];// work = arMatrixDet_MATRIX_get(ap,i, k, rowa) /arMatrixDet_MATRIX_get(ap, k, k,rowa);\r
+ for (int j = k + 1; j < dimen; j++) {\r
+ // MATRIX(ap, i, j, rowa) -= work * MATRIX(ap, k, j, rowa);\r
+ ap[i][j] -= work * ap[k][j];\r
+ }\r
+ }\r
+ }\r
+ for (int i = 0; i < dimen; i++) {\r
+ det = ap[i][i];// det *= MATRIX(ap, i, i, rowa);\r
+ }\r
+ for (int i = 0; i < is_; i++) {\r
+ det *= -1.0;\r
+ }\r
+ return det;\r
+ }\r
+\r
+ /* double arMatrixDet(ARMat *m); */\r
+ public static double arMatrixDet(NyARMat m) throws NyARException\r
+ {\r
+ NyARException.trap("動作未チェック/配列化未チェック");\r
+ if (m.row != m.clm) {\r
+ return 0.0;\r
+ }\r
+ return Det_mdet(m.getArray(), m.row, m.clm);// return mdet(m->m, m->row,m->row);\r
+ }\r
+\r
+ private final NyARVec wk_vecTridiagonalize_vec = new NyARVec(0);\r
+\r
+ private final NyARVec wk_vecTridiagonalize_vec2 = new NyARVec(0);\r
+\r
+ /**\r
+ * arVecTridiagonalize関数の代替品 a,d,e間で演算をしてる。何をどうしているかはさっぱりさっぱり\r
+ * \r
+ * @param a\r
+ * @param d\r
+ * @param e\r
+ * @param i_e_start\r
+ * 演算開始列(よくわからないけどarVecTridiagonalizeの呼び出し元でなんかしてる)\r
+ * @return\r
+ * @throws NyARException\r
+ */\r
+ private void vecTridiagonalize(NyARVec d, NyARVec e, int i_e_start)throws NyARException\r
+ {\r
+ NyARVec vec = wk_vecTridiagonalize_vec;\r
+ // double[][] a_array=a.getArray();\r
+ double s, t, p, q;\r
+ int dim;\r
+\r
+ if (this.clm != this.row) {// if(a.getClm()!=a.getRow()){\r
+ throw new NyARException();\r
+ }\r
+ if (this.clm != d.getClm()) {// if(a.getClm() != d.clm){\r
+ throw new NyARException();\r
+ }\r
+ if (this.clm != e.getClm()) {// if(a.getClm() != e.clm){\r
+ throw new NyARException();\r
+ }\r
+ dim = this.getClm();\r
+\r
+ double[] d_vec, e_vec;\r
+ d_vec = d.getArray();\r
+ e_vec = e.getArray();\r
+ double[] a_vec_k;\r
+\r
+ for (int k = 0; k < dim - 2; k++) {\r
+\r
+ a_vec_k = this.m[k];\r
+ vec.setNewArray(a_vec_k, clm);// vec=this.getRowVec(k);//double[]\r
+ // vec_array=vec.getArray();\r
+ NyARException.trap("未チェックパス");\r
+ d_vec[k] = a_vec_k[k];// d.v[k]=vec.v[k];//d.set(k,v.get(k));\r
+ // //d->v[k] = v[k];\r
+\r
+ // wv1.clm = dim-k-1;\r
+ // wv1.v = &(v[k+1]);\r
+ NyARException.trap("未チェックパス");\r
+ e_vec[k + i_e_start] = vec.vecHousehold(k + 1);// e.v[k+i_e_start]=vec.vecHousehold(k+1);//e->v[k]= arVecHousehold(&wv1);\r
+ if (e_vec[k + i_e_start] == 0.0) {// if(e.v[k+i_e_start]== 0.0){//if(e.v[k+i_e_start]== 0.0){\r
+ continue;\r
+ }\r
+\r
+ for (int i = k + 1; i < dim; i++) {\r
+ s = 0.0;\r
+ for (int j = k + 1; j < i; j++) {\r
+ NyARException.trap("未チェックのパス");\r
+ s += this.m[j][i] * a_vec_k[j];// s += a_array[j][i] *vec.v[j];//s +=a.get(j*dim+i) *v.get(j);//s +=a->m[j*dim+i] * v[j];\r
+ }\r
+ for (int j = i; j < dim; j++) {\r
+ NyARException.trap("未チェックのパス");\r
+ s += this.m[i][j] * a_vec_k[j];// s += a_array[i][j] *vec.v[j];//s +=a.get(i*dim+j) *v.get(j);//s +=a->m[i*dim+j] * v[j];\r
+ }\r
+ NyARException.trap("未チェックのパス");\r
+ d_vec[i] = s;// d.v[i]=s;//d->v[i] = s;\r
+ }\r
+\r
+ // wv1.clm = wv2.clm = dim-k-1;\r
+ // wv1.v = &(v[k+1]);\r
+ // wv2.v = &(d->v[k+1]);\r
+ a_vec_k = this.m[k];\r
+ vec.setNewArray(a_vec_k, clm);// vec=this.getRowVec(k);\r
+ // vec_array=vec.getArray();\r
+ NyARException.trap("未チェックパス");\r
+ t = vec.vecInnerproduct(d, k + 1) / 2;\r
+ for (int i = dim - 1; i > k; i--) {\r
+ NyARException.trap("未チェックパス");\r
+ p = a_vec_k[i];// p = v.get(i);//p = v[i];\r
+ d_vec[i] -= t * p;\r
+ q = d_vec[i];// d.v[i]-=t*p;q=d.v[i];//q = d->v[i] -= t*p;\r
+ for (int j = i; j < dim; j++) {\r
+ NyARException.trap("未チェックパス");\r
+ this.m[i][j] -= p * (d_vec[j] + q * a_vec_k[j]);// a.m[i][j]-=p*(d.v[j] +q*vec.v[j]);//a->m[i*dim+j] -=p*(d->v[j]) + q*v[j];\r
+ }\r
+ }\r
+ }\r
+\r
+ if (dim >= 2) {\r
+ d_vec[dim - 2] = this.m[dim - 2][dim - 2];// d.v[dim-2]=a.m[dim-2][dim-2];//d->v[dim-2]=a->m[(dim-2)*dim+(dim-2)];\r
+ e_vec[dim - 2 + i_e_start] = this.m[dim - 2][dim - 1];// e.v[dim-2+i_e_start]=a.m[dim-2][dim-1];//e->v[dim-2] = a->m[(dim-2)*dim+(dim-1)];\r
+ }\r
+\r
+ if (dim >= 1) {\r
+ d_vec[dim - 1] = this.m[dim - 1][dim - 1];// d.v[dim-1]=a_array[dim-1][dim-1];//d->v[dim-1] =a->m[(dim-1)*dim+(dim-1)];\r
+ }\r
+ NyARVec vec2 = this.wk_vecTridiagonalize_vec2;\r
+ for (int k = dim - 1; k >= 0; k--) {\r
+ a_vec_k = this.m[k];\r
+ vec.setNewArray(a_vec_k, clm);// vec=this.getRowVec(k);//v =a.getPointer(k*dim);//v = &(a->m[k*dim]);\r
+ if (k < dim - 2) {\r
+ for (int i = k + 1; i < dim; i++) {\r
+ // wv1.clm = wv2.clm = dim-k-1;\r
+ // wv1.v = &(v[k+1]);\r
+ // wv2.v = &(a->m[i*dim+k+1]);\r
+ vec2.setNewArray(this.m[i], clm);// vec2=this.getRowVec(i);\r
+\r
+ t = vec.vecInnerproduct(vec2, k + 1);\r
+ for (int j = k + 1; j < dim; j++) {\r
+ NyARException.trap("未チェックパス");\r
+ this.m[i][j] -= t * a_vec_k[j];// a_array[i][j]-=t*vec.v[j];//a.subValue(i*dim+j,t*v.get(j));//a->m[i*dim+j]-= t * v[j];\r
+ }\r
+ }\r
+ }\r
+ for (int i = 0; i < dim; i++) {\r
+ a_vec_k[i] = 0.0;// v.set(i,0.0);//v[i] = 0.0;\r
+ }\r
+ a_vec_k[k] = 1;// v.set(k,1);//v[k] = 1;\r
+ }\r
+ return;\r
+ }\r
+}
\ No newline at end of file