tags/2.3.1/src/jp/nyatla/nyartoolkit/core/NyARMat.java

   1 /* \r
   2  * PROJECT: NyARToolkit\r
   3  * --------------------------------------------------------------------------------\r
   4  * This work is based on the original ARToolKit developed by\r
   5  *   Hirokazu Kato\r
   6  *   Mark Billinghurst\r
   7  *   HITLab, University of Washington, Seattle\r
   8  * http://www.hitl.washington.edu/artoolkit/\r
   9  *\r
  10  * The NyARToolkit is Java version ARToolkit class library.\r
  11  * Copyright (C)2008 R.Iizuka\r
  12  *\r
  13  * This program is free software; you can redistribute it and/or\r
  14  * modify it under the terms of the GNU General Public License\r
  15  * as published by the Free Software Foundation; either version 2\r
  16  * of the License, or (at your option) any later version.\r
  17  * \r
  18  * This program is distributed in the hope that it will be useful,\r
  19  * but WITHOUT ANY WARRANTY; without even the implied warranty of\r
  20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\r
  21  * GNU General Public License for more details.\r
  22  * \r
  23  * You should have received a copy of the GNU General Public License\r
  24  * along with this framework; if not, write to the Free Software\r
  25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA\r
  26  * \r
  27  * For further information please contact.\r
  28  *      http://nyatla.jp/nyatoolkit/\r
  29  *      <airmail(at)ebony.plala.or.jp>\r
  30  * \r
  31  */\r
  32 package jp.nyatla.nyartoolkit.core;\r
  33 \r
  34 import jp.nyatla.nyartoolkit.NyARException;\r
  35 \r
  36 \r
  37 \r
  38 /**\r
  39  * ARMat構造体に対応するクラス typedef struct { double *m; int row; int clm; }ARMat;\r
  40  * \r
  41  */\r
  42 public class NyARMat\r
  43 {\r
  44         /**\r
  45          * 配列サイズと行列サイズは必ずしも一致しないことに注意 返された配列のサイズを行列の大きさとして使わないこと！\r
  46          * \r
  47          */\r
  48         protected double[][] _m;\r
  49         private int[] __matrixSelfInv_nos;\r
  50 \r
  51         protected int clm;\r
  52         protected int row;\r
  53 \r
  54         /**\r
  55          * デフォルトコンストラクタは機能しません。\r
  56          * \r
  57          * @throws NyARException\r
  58          */\r
  59         protected NyARMat() throws NyARException\r
  60         {\r
  61                 throw new NyARException();\r
  62         }\r
  63 \r
  64         public NyARMat(int i_row, int i_clm)\r
  65         {\r
  66                 this._m = new double[i_row][i_clm];\r
  67                 this.clm = i_clm;\r
  68                 this.row = i_row;\r
  69                 this.__matrixSelfInv_nos=new int[i_row];\r
  70                 return;\r
  71         }\r
  72         \r
  73         /**\r
  74          * 行列の列数を返します。\r
  75          * @return\r
  76          */\r
  77         public int getClm()\r
  78         {\r
  79                 return clm;\r
  80         }\r
  81         /**\r
  82          * 行列の行数を返します。\r
  83          * @return\r
  84          */\r
  85         public int getRow()\r
  86         {\r
  87                 return row;\r
  88         }       \r
  89         \r
  90         /**\r
  91          * 行列のサイズを変更します。\r
  92          * 実行後、行列の各値は不定になります。\r
  93          * @param i_row\r
  94          * @param i_clm\r
  95          */\r
  96         public void realloc(int i_row, int i_clm)\r
  97         {\r
  98                 if (i_row <= this._m.length && i_clm <= this._m[0].length) {\r
  99                         // 十分な配列があれば何もしない。\r
 100                 } else {\r
 101                         // 不十分なら取り直す。\r
 102                         this._m = new double[i_row][i_clm];\r
 103                         this.__matrixSelfInv_nos=new int[i_row];\r
 104                 }\r
 105                 this.clm = i_clm;\r
 106                 this.row = i_row;\r
 107                 return;\r
 108         }\r
 109         \r
 110         \r
 111         \r
 112         /**\r
 113          * i_mat_aとi_mat_bの積を計算して、thisへ格納します。\r
 114          * @param i_mat_a\r
 115          * @param i_mat_b\r
 116          * @throws NyARException\r
 117          */\r
 118         public void matrixMul(NyARMat i_mat_a, NyARMat i_mat_b) throws NyARException\r
 119         {\r
 120                 assert i_mat_a.clm == i_mat_b.row && this.row==i_mat_a.row && this.clm==i_mat_b.clm;\r
 121 \r
 122                 double w;\r
 123                 int r, c, i;\r
 124                 double[][] am = i_mat_a._m, bm = i_mat_b._m, dm = this._m;\r
 125                 // For順変更禁止\r
 126                 for (r = 0; r < this.row; r++) {\r
 127                         for (c = 0; c < this.clm; c++) {\r
 128                                 w = 0.0;\r
 129                                 for (i = 0; i < i_mat_a.clm; i++) {\r
 130                                         w += am[r][i] * bm[i][c];\r
 131                                 }\r
 132                                 dm[r][c] = w;\r
 133                         }\r
 134                 }\r
 135                 return;\r
 136         }       \r
 137         \r
 138         /**\r
 139          * 逆行列を計算して、thisへ格納します。\r
 140          * @throws NyARException\r
 141          */\r
 142         public boolean matrixSelfInv() throws NyARException\r
 143         {\r
 144                 double[][] ap = this._m;\r
 145                 int dimen = this.row;\r
 146                 int dimen_1 = dimen - 1;\r
 147                 double[] ap_n, ap_ip, ap_i;// wap;\r
 148                 int j, ip, nwork;\r
 149                 int[] nos = __matrixSelfInv_nos;//ワーク変数\r
 150                 // double epsl;\r
 151                 double p, pbuf, work;\r
 152 \r
 153                 /* check size */\r
 154                 switch (dimen) {\r
 155                 case 0:\r
 156                         throw new NyARException();\r
 157                 case 1:\r
 158                         ap[0][0] = 1.0 / ap[0][0];// *ap = 1.0 / (*ap);\r
 159                         return true;/* 1 dimension */\r
 160                 }\r
 161 \r
 162                 for (int n = 0; n < dimen; n++) {\r
 163                         nos[n] = n;\r
 164                 }\r
 165 \r
 166                 /*\r
 167                  * nyatla memo ipが定まらないで計算が行われる場合があるので挿入。 ループ内で0初期化していいかが判らない。\r
 168                  */\r
 169                 ip = 0;\r
 170                 // For順変更禁止\r
 171                 for (int n = 0; n < dimen; n++) {\r
 172                         ap_n = ap[n];// wcp = ap + n * rowa;\r
 173                         p = 0.0;\r
 174                         for (int i = n; i < dimen; i++) {\r
 175                                 if (p < (pbuf = Math.abs(ap[i][0]))) {\r
 176                                         p = pbuf;\r
 177                                         ip = i;\r
 178                                 }\r
 179                         }\r
 180                         // if (p <= matrixSelfInv_epsl){\r
 181                         if (p == 0.0) {\r
 182                                 return false;\r
 183                                 // throw new NyARException();\r
 184                         }\r
 185 \r
 186                         nwork = nos[ip];\r
 187                         nos[ip] = nos[n];\r
 188                         nos[n] = nwork;\r
 189 \r
 190                         ap_ip = ap[ip];\r
 191                         for (j = 0; j < dimen; j++) {// for(j = 0, wap = ap + ip * rowa,\r
 192                                                                                         // wbp = wcp; j < dimen ; j++) {\r
 193                                 work = ap_ip[j]; // work = *wap;\r
 194                                 ap_ip[j] = ap_n[j];\r
 195                                 ap_n[j] = work;\r
 196                         }\r
 197 \r
 198                         work = ap_n[0];\r
 199                         for (j = 0; j < dimen_1; j++) {\r
 200                                 ap_n[j] = ap_n[j + 1] / work;// *wap = *(wap + 1) / work;\r
 201                         }\r
 202                         ap_n[j] = 1.0 / work;// *wap = 1.0 / work;\r
 203                         for (int i = 0; i < dimen; i++) {\r
 204                                 if (i != n) {\r
 205                                         ap_i = ap[i];// wap = ap + i * rowa;\r
 206                                         work = ap_i[0];\r
 207                                         for (j = 0; j < dimen_1; j++) {// for(j = 1, wbp = wcp,work = *wap;j < dimen ;j++, wap++, wbp++)\r
 208                                                 ap_i[j] = ap_i[j + 1] - work * ap_n[j];// wap = *(wap +1) - work *(*wbp);\r
 209                                         }\r
 210                                         ap_i[j] = -work * ap_n[j];// *wap = -work * (*wbp);\r
 211                                 }\r
 212                         }\r
 213                 }\r
 214 \r
 215                 for (int n = 0; n < dimen; n++) {\r
 216                         for (j = n; j < dimen; j++) {\r
 217                                 if (nos[j] == n) {\r
 218                                         break;\r
 219                                 }\r
 220                         }\r
 221                         nos[j] = nos[n];\r
 222                         for (int i = 0; i < dimen; i++) {\r
 223                                 ap_i = ap[i];\r
 224                                 work = ap_i[j];// work = *wap;\r
 225                                 ap_i[j] = ap_i[n];// *wap = *wbp;\r
 226                                 ap_i[n] = work;// *wbp = work;\r
 227                         }\r
 228                 }\r
 229                 return true;\r
 230         }\r
 231 \r
 232         /**\r
 233          * 行列をゼロクリアする。\r
 234          */\r
 235         public void zeroClear()\r
 236         {\r
 237                 int i, i2;\r
 238                 // For順変更OK\r
 239                 for (i = row - 1; i >= 0; i--) {\r
 240                         for (i2 = clm - 1; i2 >= 0; i2--) {\r
 241                                 _m[i][i2] = 0.0;\r
 242                         }\r
 243                 }\r
 244         }\r
 245 \r
 246         /**\r
 247          * i_copy_fromの内容を自分自身にコピーします。 高さ・幅は同一で無いと失敗します。\r
 248          * \r
 249          * @param i_copy_from\r
 250          */\r
 251         public void copyFrom(NyARMat i_copy_from) throws NyARException\r
 252         {\r
 253                 // サイズ確認\r
 254                 if (this.row != i_copy_from.row || this.clm != i_copy_from.clm) {\r
 255                         throw new NyARException();\r
 256                 }\r
 257                 // 値コピー\r
 258                 for (int r = this.row - 1; r >= 0; r--) {\r
 259                         for (int c = this.clm - 1; c >= 0; c--) {\r
 260                                 this._m[r][c] = i_copy_from._m[r][c];\r
 261                         }\r
 262                 }\r
 263         }\r
 264 \r
 265         public double[][] getArray()\r
 266         {\r
 267                 return _m;\r
 268         }\r
 269 \r
 270 \r
 271         /**\r
 272          * sourceの転置行列をdestに得る。arMatrixTrans()の代替品\r
 273          * \r
 274          * @param dest\r
 275          * @param source\r
 276          * @return\r
 277          */\r
 278         public static void matrixTrans(NyARMat dest, NyARMat source) throws NyARException\r
 279         {\r
 280                 if (dest.row != source.clm || dest.clm != source.row) {\r
 281                         throw new NyARException();\r
 282                 }\r
 283                 NyARException.trap("未チェックのパス");\r
 284                 // For順変更禁止\r
 285                 for (int r = 0; r < dest.row; r++) {\r
 286                         for (int c = 0; c < dest.clm; c++) {\r
 287                                 dest._m[r][c] = source._m[c][r];\r
 288                         }\r
 289                 }\r
 290         }\r
 291 \r
 292         /**\r
 293          * unitを単位行列に初期化する。arMatrixUnitの代替品\r
 294          * \r
 295          * @param unit\r
 296          */\r
 297         public static void matrixUnit(NyARMat unit) throws NyARException\r
 298         {\r
 299                 if (unit.row != unit.clm) {\r
 300                         throw new NyARException();\r
 301                 }\r
 302                 NyARException.trap("未チェックのパス");\r
 303                 // For順変更禁止\r
 304                 for (int r = 0; r < unit.getRow(); r++) {\r
 305                         for (int c = 0; c < unit.getClm(); c++) {\r
 306                                 if (r == c) {\r
 307                                         unit._m[r][c] = 1.0;\r
 308                                 } else {\r
 309                                         unit._m[r][c] = 0.0;\r
 310                                 }\r
 311                         }\r
 312                 }\r
 313         }\r
 314 \r
 315         /**\r
 316          * sourceの内容を自身に複製する。 Optimized 2008.04.19\r
 317          * \r
 318          * @param i_source\r
 319          * @return\r
 320          */\r
 321         public void matrixDup(NyARMat i_source) throws NyARException\r
 322         {\r
 323                 // 自身の配列サイズを相手のそれより大きいことを保障する。\r
 324                 this.realloc(i_source.row, i_source.clm);\r
 325                 // 内容を転写\r
 326                 int r, c;\r
 327                 double[][] src_m, dest_m;\r
 328                 src_m = i_source._m;\r
 329                 dest_m = this._m;\r
 330                 // コピーはFor順を変えてもOK\r
 331                 for (r = this.row - 1; r >= 0; r--) {\r
 332                         for (c = this.clm - 1; c >= 0; c--) {\r
 333                                 dest_m[r][c] = src_m[r][c];\r
 334                         }\r
 335                 }\r
 336         }\r
 337 \r
 338         public NyARMat matrixAllocDup() throws NyARException\r
 339         {\r
 340                 NyARMat result = new NyARMat(this.row, this.clm);\r
 341                 // コピー\r
 342                 int r, c;\r
 343                 double[][] dest_m, src_m;\r
 344                 dest_m = result._m;\r
 345                 src_m = this._m;\r
 346                 // コピーはFor順を変えてもOK\r
 347                 for (r = this.row - 1; r >= 0; r--) {\r
 348                         for (c = this.clm - 1; c >= 0; c--) {\r
 349                                 dest_m[r][c] = src_m[r][c];\r
 350                         }\r
 351                 }\r
 352                 return result;\r
 353         }\r
 354 \r
 355 \r
 356         private static final double PCA_EPS = 1e-6; // #define EPS 1e-6\r
 357 \r
 358         private static final int PCA_MAX_ITER = 100; // #define MAX_ITER 100\r
 359 \r
 360         private static final double PCA_VZERO = 1e-16; // #define VZERO 1e-16\r
 361 \r
 362         /**\r
 363          * static int EX( ARMat *input, ARVec *mean )の代替関数 Optimize:STEP:[144->110]\r
 364          * \r
 365          * @param input\r
 366          * @param mean\r
 367          * @return\r
 368          * @throws NyARException\r
 369          */\r
 370         private void PCA_EX(NyARVec mean) throws NyARException\r
 371         {\r
 372                 int lrow, lclm;\r
 373                 int i, i2;\r
 374                 lrow = this.row;\r
 375                 lclm = this.clm;\r
 376                 double[][] lm = this._m;\r
 377 \r
 378                 if (lrow <= 0 || lclm <= 0) {\r
 379                         throw new NyARException();\r
 380                 }\r
 381                 if (mean.getClm() != lclm) {\r
 382                         throw new NyARException();\r
 383                 }\r
 384                 // double[] mean_array=mean.getArray();\r
 385                 // mean.zeroClear();\r
 386                 final double[] mean_array = mean.getArray();\r
 387                 double w;\r
 388                 // For順変更禁止\r
 389                 for (i2 = 0; i2 < lclm; i2++) {\r
 390                         w = 0.0;\r
 391                         for (i = 0; i < lrow; i++) {\r
 392                                 // *(v++) += *(m++);\r
 393                                 w += lm[i][i2];\r
 394                         }\r
 395                         mean_array[i2] = w / lrow;// mean->v[i] /= row;\r
 396                 }\r
 397         }\r
 398 \r
 399         /**\r
 400          * static int CENTER( ARMat *inout, ARVec *mean )の代替関数\r
 401          * \r
 402          * @param inout\r
 403          * @param mean\r
 404          * @return\r
 405          */\r
 406         private static void PCA_CENTER(NyARMat inout, NyARVec mean)throws NyARException\r
 407         {\r
 408                 double[] v;\r
 409                 int row, clm;\r
 410 \r
 411                 row = inout.row;\r
 412                 clm = inout.clm;\r
 413                 if (mean.getClm() != clm) {\r
 414                         throw new NyARException();\r
 415                 }\r
 416                 double[][] im = inout._m;\r
 417                 double[] im_i;\r
 418                 double w0, w1;\r
 419                 v = mean.getArray();\r
 420                 // 特にパフォーマンスが劣化するclm=1と2ときだけ、別パスで処理します。\r
 421                 switch (clm) {\r
 422                 case 1:\r
 423                         w0 = v[0];\r
 424                         for (int i = 0; i < row; i++) {\r
 425                                 im[i][0] -= w0;\r
 426                         }\r
 427                         break;\r
 428                 case 2:\r
 429                         w0 = v[0];\r
 430                         w1 = v[1];\r
 431                         for (int i = 0; i < row; i++) {\r
 432                                 im_i = im[i];\r
 433                                 im_i[0] -= w0;\r
 434                                 im_i[1] -= w1;\r
 435                         }\r
 436                         break;\r
 437                 default:\r
 438                         for (int i = 0; i < row; i++) {\r
 439                                 im_i = im[i];\r
 440                                 for (int j = 0; j < clm; j++) {\r
 441                                         // *(m++) -= *(v++);\r
 442                                         im_i[j] -= v[j];\r
 443                                 }\r
 444                         }\r
 445                         break;\r
 446                 }\r
 447                 return;\r
 448         }\r
 449 \r
 450         /**\r
 451          * int x_by_xt( ARMat *input, ARMat *output )の代替関数\r
 452          * \r
 453          * @param input\r
 454          * @param output\r
 455          * @throws NyARException\r
 456          */\r
 457         private static void PCA_x_by_xt(NyARMat input, NyARMat output) throws NyARException\r
 458         {\r
 459                 NyARException.trap("動作未チェック/配列化未チェック");\r
 460                 int row, clm;\r
 461                 // double[][] out;\r
 462                 double[] in1, in2;\r
 463 \r
 464                 NyARException.trap("未チェックのパス");\r
 465                 row = input.row;\r
 466                 clm = input.clm;\r
 467                 NyARException.trap("未チェックのパス");\r
 468                 if (output.row != row || output.clm != row) {\r
 469                         throw new NyARException();\r
 470                 }\r
 471 \r
 472                 // out = output.getArray();\r
 473                 for (int i = 0; i < row; i++) {\r
 474                         for (int j = 0; j < row; j++) {\r
 475                                 if (j < i) {\r
 476                                         NyARException.trap("未チェックのパス");\r
 477                                         output._m[i][j] = output._m[j][i];// *out =\r
 478                                                                                                         // output->m[j*row+i];\r
 479                                 } else {\r
 480                                         NyARException.trap("未チェックのパス");\r
 481                                         in1 = input._m[i];// input.getRowArray(i);//in1 = &(input->m[clm*i]);\r
 482                                         in2 = input._m[j];// input.getRowArray(j);//in2 = &(input->m[clm*j]);\r
 483                                         output._m[i][j] = 0;// *out = 0.0;\r
 484                                         for (int k = 0; k < clm; k++) {\r
 485                                                 output._m[i][j] += (in1[k] * in2[k]);// *out += *(in1++)\r
 486                                                                                                                         // * *(in2++);\r
 487                                         }\r
 488                                 }\r
 489                                 // out.incPtr();\r
 490                         }\r
 491                 }\r
 492         }\r
 493 \r
 494         /**\r
 495          * static int xt_by_x( ARMat *input, ARMat *output )の代替関数\r
 496          * Optimize:2008.04.19\r
 497          * \r
 498          * @param input\r
 499          * @param i_output\r
 500          * @throws NyARException\r
 501          */\r
 502         private static void PCA_xt_by_x(NyARMat input, NyARMat i_output) throws NyARException\r
 503         {\r
 504                 double[] in_;\r
 505                 int row, clm;\r
 506 \r
 507                 row = input.row;\r
 508                 clm = input.clm;\r
 509                 if (i_output.row != clm || i_output.clm != clm) {\r
 510                         throw new NyARException();\r
 511                 }\r
 512 \r
 513                 int k, j;\r
 514                 double[][] out_m = i_output._m;\r
 515                 double w;\r
 516                 for (int i = 0; i < clm; i++) {\r
 517                         for (j = 0; j < clm; j++) {\r
 518                                 if (j < i) {\r
 519                                         out_m[i][j] = out_m[j][i];// *out = output->m[j*clm+i];\r
 520                                 } else {\r
 521                                         w = 0.0;// *out = 0.0;\r
 522                                         for (k = 0; k < row; k++) {\r
 523                                                 in_ = input._m[k];// in=input.getRowArray(k);\r
 524                                                 w += (in_[i] * in_[j]);// *out += *in1 * *in2;\r
 525                                         }\r
 526                                         out_m[i][j] = w;\r
 527                                 }\r
 528                         }\r
 529                 }\r
 530         }\r
 531 \r
 532         private final NyARVec wk_PCA_QRM_ev = new NyARVec(1);\r
 533 \r
 534         /**\r
 535          * static int QRM( ARMat *a, ARVec *dv )の代替関数\r
 536          * \r
 537          * @param a\r
 538          * @param dv\r
 539          * @throws NyARException\r
 540          */\r
 541         private void PCA_QRM(NyARVec dv) throws NyARException\r
 542         {\r
 543                 double w, t, s, x, y, c;\r
 544                 int dim, iter;\r
 545                 double[] dv_array = dv.getArray();\r
 546 \r
 547                 dim = this.row;\r
 548                 if (dim != this.clm || dim < 2) {\r
 549                         throw new NyARException();\r
 550                 }\r
 551                 if (dv.getClm() != dim) {\r
 552                         throw new NyARException();\r
 553                 }\r
 554 \r
 555                 NyARVec ev = this.wk_PCA_QRM_ev;\r
 556                 ev.realloc(dim);\r
 557                 double[] ev_array = ev.getArray();\r
 558                 if (ev == null) {\r
 559                         throw new NyARException();\r
 560                 }\r
 561                 final double[][] L_m = this._m;\r
 562                 this.vecTridiagonalize(dv, ev, 1);\r
 563 \r
 564                 ev_array[0] = 0.0;// ev->v[0] = 0.0;\r
 565                 for (int h = dim - 1; h > 0; h--) {\r
 566                         int j = h;\r
 567                         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
 568                                 // j--;\r
 569                                 j--;\r
 570                         }\r
 571                         if (j == h) {\r
 572                                 continue;\r
 573                         }\r
 574                         iter = 0;\r
 575                         do {\r
 576                                 iter++;\r
 577                                 if (iter > PCA_MAX_ITER) {\r
 578                                         break;\r
 579                                 }\r
 580                                 w = (dv_array[h - 1] - dv_array[h]) / 2;// w = (dv->v[h-1] -dv->v[h]) / 2;//ここ？\r
 581                                 t = ev_array[h] * ev_array[h];// t = ev->v[h] * ev->v[h];\r
 582                                 s = Math.sqrt(w * w + t);\r
 583                                 if (w < 0) {\r
 584                                         s = -s;\r
 585                                 }\r
 586                                 x = dv_array[j] - dv_array[h] + t / (w + s);// x = dv->v[j] -dv->v[h] +t/(w+s);\r
 587                                 y = ev_array[j + 1];// y = ev->v[j+1];\r
 588                                 for (int k = j; k < h; k++) {\r
 589                                         if (Math.abs(x) >= Math.abs(y)) {\r
 590                                                 if (Math.abs(x) > PCA_VZERO) {\r
 591                                                         t = -y / x;\r
 592                                                         c = 1 / Math.sqrt(t * t + 1);\r
 593                                                         s = t * c;\r
 594                                                 } else {\r
 595                                                         c = 1.0;\r
 596                                                         s = 0.0;\r
 597                                                 }\r
 598                                         } else {\r
 599                                                 t = -x / y;\r
 600                                                 s = 1.0 / Math.sqrt(t * t + 1);\r
 601                                                 c = t * s;\r
 602                                         }\r
 603                                         w = dv_array[k] - dv_array[k + 1];// w = dv->v[k] -dv->v[k+1];\r
 604                                         t = (w * s + 2 * c * ev_array[k + 1]) * s;// t = (w * s +2 * c *ev->v[k+1]) *s;\r
 605                                         dv_array[k] -= t;// dv->v[k] -= t;\r
 606                                         dv_array[k + 1] += t;// dv->v[k+1] += t;\r
 607                                         if (k > j) {\r
 608                                                 NyARException.trap("未チェックパス");\r
 609                                                 {\r
 610                                                         ev_array[k] = c * ev_array[k] - s * y;// ev->v[k]= c *ev->v[k]- s * y;\r
 611                                                 }\r
 612                                         }\r
 613                                         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
 614 \r
 615                                         for (int i = 0; i < dim; i++) {\r
 616                                                 x = L_m[k][i];// x = a->m[k*dim+i];\r
 617                                                 y = L_m[k + 1][i];// y = a->m[(k+1)*dim+i];\r
 618                                                 L_m[k][i] = c * x - s * y;// a->m[k*dim+i] = c * x - s* y;\r
 619                                                 L_m[k + 1][i] = s * x + c * y;// a->m[(k+1)*dim+i] = s* x + c * y;\r
 620                                         }\r
 621                                         if (k < h - 1) {\r
 622                                                 NyARException.trap("未チェックパス");\r
 623                                                 {\r
 624                                                         x = ev_array[k + 1];// x = ev->v[k+1];\r
 625                                                         y = -s * ev_array[k + 2];// y = -s * ev->v[k+2];\r
 626                                                         ev_array[k + 2] *= c;// ev->v[k+2] *= c;\r
 627                                                 }\r
 628                                         }\r
 629                                 }\r
 630                         } while (Math.abs(ev_array[h]) > PCA_EPS\r
 631                                         * (Math.abs(dv_array[h - 1]) + Math.abs(dv_array[h])));\r
 632                 }\r
 633                 for (int k = 0; k < dim - 1; k++) {\r
 634                         int h = k;\r
 635                         t = dv_array[h];// t = dv->v[h];\r
 636                         for (int i = k + 1; i < dim; i++) {\r
 637                                 if (dv_array[i] > t) {// if( dv->v[i] > t ) {\r
 638                                         h = i;\r
 639                                         t = dv_array[h];// t = dv->v[h];\r
 640                                 }\r
 641                         }\r
 642                         dv_array[h] = dv_array[k];// dv->v[h] = dv->v[k];\r
 643                         dv_array[k] = t;// dv->v[k] = t;\r
 644                         this.flipRow(h, k);\r
 645                 }\r
 646         }\r
 647 \r
 648         /**\r
 649          * i_row_1番目の行と、i_row_2番目の行を入れ替える。\r
 650          * \r
 651          * @param i_row_1\r
 652          * @param i_row_2\r
 653          */\r
 654         private void flipRow(int i_row_1, int i_row_2)\r
 655         {\r
 656                 int i;\r
 657                 double w;\r
 658                 double[] r1 = this._m[i_row_1], r2 = this._m[i_row_2];\r
 659                 // For順変更OK\r
 660                 for (i = clm - 1; i >= 0; i--) {\r
 661                         w = r1[i];\r
 662                         r1[i] = r2[i];\r
 663                         r2[i] = w;\r
 664                 }\r
 665         }\r
 666 \r
 667         /**\r
 668          * static int EV_create( ARMat *input, ARMat *u, ARMat *output, ARVec *ev\r
 669          * )の代替関数\r
 670          * \r
 671          * @param input\r
 672          * @param u\r
 673          * @param output\r
 674          * @param ev\r
 675          * @throws NyARException\r
 676          */\r
 677         private static void PCA_EV_create(NyARMat input, NyARMat u, NyARMat output,NyARVec ev) throws NyARException\r
 678         {\r
 679                 NyARException.trap("未チェックのパス");\r
 680                 int row, clm;\r
 681                 row = input.row;// row = input->row;\r
 682                 clm = input.clm;// clm = input->clm;\r
 683                 if (row <= 0 || clm <= 0) {\r
 684                         throw new NyARException();\r
 685                 }\r
 686                 if (u.row != row || u.clm != row) {// if( u->row != row || u->clm !=\r
 687                                                                                         // row ){\r
 688                         throw new NyARException();\r
 689                 }\r
 690                 if (output.row != row || output.clm != clm) {// if( output->row !=\r
 691                                                                                                                 // row || output->clm !=\r
 692                                                                                                                 // clm ){\r
 693                         throw new NyARException();\r
 694                 }\r
 695                 if (ev.getClm() != row) {// if( ev->clm != row ){\r
 696                         throw new NyARException();\r
 697                 }\r
 698                 double[][] m, in_;\r
 699                 double[] m1, ev_array;\r
 700                 double sum, work;\r
 701 \r
 702                 NyARException.trap("未チェックのパス");\r
 703                 m = output._m;// m = output->m;\r
 704                 in_ = input._m;\r
 705                 int i;\r
 706                 ev_array = ev.getArray();\r
 707                 for (i = 0; i < row; i++) {\r
 708                         NyARException.trap("未チェックのパス");\r
 709                         if (ev_array[i] < PCA_VZERO) {// if( ev->v[i] < VZERO ){\r
 710                                 break;\r
 711                         }\r
 712                         NyARException.trap("未チェックのパス");\r
 713                         work = 1 / Math.sqrt(Math.abs(ev_array[i]));// work = 1 /\r
 714                                                                                                                 // sqrt(fabs(ev->v[i]));\r
 715                         for (int j = 0; j < clm; j++) {\r
 716                                 sum = 0.0;\r
 717                                 m1 = u._m[i];// m1 = &(u->m[i*row]);\r
 718                                 // m2=input.getPointer(j);//m2 = &(input->m[j]);\r
 719                                 for (int k = 0; k < row; k++) {\r
 720                                         sum += m1[k] + in_[k][j];// sum += *m1 * *m2;\r
 721                                         // m1.incPtr(); //m1++;\r
 722                                         // m2.addPtr(clm);//m2 += clm;\r
 723                                 }\r
 724                                 m1[j] = sum * work;// *(m++) = sum * work;\r
 725                                 // {//*(m++) = sum * work;\r
 726                                 // m.set(sum * work);\r
 727                                 // m.incPtr();}\r
 728                         }\r
 729                 }\r
 730                 for (; i < row; i++) {\r
 731                         NyARException.trap("未チェックのパス");\r
 732                         ev_array[i] = 0.0;// ev->v[i] = 0.0;\r
 733                         for (int j = 0; j < clm; j++) {\r
 734                                 m[i][j] = 0.0;\r
 735                                 // m.set(0.0);//*(m++) = 0.0;\r
 736                                 // m.incPtr();\r
 737                         }\r
 738                 }\r
 739         }\r
 740 \r
 741         private NyARMat wk_PCA_PCA_u = null;\r
 742 \r
 743         /**\r
 744          * static int PCA( ARMat *input, ARMat *output, ARVec *ev )\r
 745          * \r
 746          * @param output\r
 747          * @param o_ev\r
 748          * @throws NyARException\r
 749          */\r
 750         private void PCA_PCA(NyARMat o_output, NyARVec o_ev) throws NyARException\r
 751         {\r
 752 \r
 753                 int l_row, l_clm, min;\r
 754                 double[] ev_array = o_ev.getArray();\r
 755 \r
 756                 l_row = this.row;// row = input->row;\r
 757                 l_clm = this.clm;// clm = input->clm;\r
 758                 min = (l_clm < l_row) ? l_clm : l_row;\r
 759                 if (l_row < 2 || l_clm < 2) {\r
 760                         throw new NyARException();\r
 761                 }\r
 762                 if (o_output.clm != this.clm) {// if( output->clm != input->clm ){\r
 763                         throw new NyARException();\r
 764                 }\r
 765                 if (o_output.row != min) {// if( output->row != min ){\r
 766                         throw new NyARException();\r
 767                 }\r
 768                 if (o_ev.getClm() != min) {// if( ev->clm != min ){\r
 769                         throw new NyARException();\r
 770                 }\r
 771 \r
 772                 NyARMat u;// u =new NyARMat( min, min );\r
 773                 if (this.wk_PCA_PCA_u == null) {\r
 774                         u = new NyARMat(min, min);\r
 775                         this.wk_PCA_PCA_u = u;\r
 776                 } else {\r
 777                         u = this.wk_PCA_PCA_u;\r
 778                         u.realloc(min, min);\r
 779                 }\r
 780 \r
 781                 if (l_row < l_clm) {\r
 782                         NyARException.trap("未チェックのパス");\r
 783                         PCA_x_by_xt(this, u);// if(x_by_xt( input, u ) < 0 ) {\r
 784                 } else {\r
 785                         PCA_xt_by_x(this, u);// if(xt_by_x( input, u ) < 0 ) {\r
 786                 }\r
 787                 u.PCA_QRM(o_ev);\r
 788 \r
 789                 double[][] m1, m2;\r
 790                 if (l_row < l_clm) {\r
 791                         NyARException.trap("未チェックのパス");\r
 792                         PCA_EV_create(this, u, o_output, o_ev);\r
 793                 } else {\r
 794                         m1 = u._m;// m1 = u->m;\r
 795                         m2 = o_output._m;// m2 = output->m;\r
 796                         int i;\r
 797                         for (i = 0; i < min; i++) {\r
 798                                 if (ev_array[i] < PCA_VZERO) {// if( ev->v[i] < VZERO ){\r
 799                                         break;\r
 800                                 }\r
 801                                 for (int j = 0; j < min; j++) {\r
 802                                         m2[i][j] = m1[i][j];// *(m2++) = *(m1++);\r
 803                                 }\r
 804                         }\r
 805                         for (; i < min; i++) {// for( ; i < min; i++){\r
 806                                 // コードを見た限りあってそうだからコメントアウト(2008/03/26)NyARException.trap("未チェックのパス");\r
 807                                 ev_array[i] = 0.0;// ev->v[i] = 0.0;\r
 808                                 for (int j = 0; j < min; j++) {\r
 809                                         m2[i][j] = 0.0;// *(m2++) = 0.0;\r
 810                                 }\r
 811                         }\r
 812                 }\r
 813         }\r
 814         /**\r
 815          * 主成分分析を実行して、結果をthisと引数へ格納します。\r
 816          * @param o_evec\r
 817          * @param o_ev\r
 818          * @param o_mean\r
 819          * @throws NyARException\r
 820          */\r
 821         public void pca(NyARMat o_evec, NyARVec o_ev, NyARVec o_mean)throws NyARException\r
 822         {\r
 823                 final double l_row = this.row;// row = input->row;\r
 824                 final double l_clm = this.clm;// clm = input->clm;\r
 825                 final double check=(l_row < l_clm) ? l_row : l_clm;\r
 826                 \r
 827                 assert l_row >= 2 || l_clm >= 2;\r
 828                 assert o_evec.clm == l_clm && o_evec.row == check;\r
 829                 assert o_ev.getClm() == check;\r
 830                 assert o_mean.getClm() == l_clm;\r
 831                 \r
 832                 final double srow = Math.sqrt((double) l_row);\r
 833                 PCA_EX(o_mean);\r
 834 \r
 835                 PCA_CENTER(this, o_mean);\r
 836 \r
 837                 int i, j;\r
 838                 // For順変更OK\r
 839                 for (i = 0; i < l_row; i++) {\r
 840                         for (j = 0; j < l_clm; j++) {\r
 841                                 this._m[i][j] /= srow;// work->m[i] /= srow;\r
 842                         }\r
 843                 }\r
 844 \r
 845                 PCA_PCA(o_evec, o_ev);\r
 846 \r
 847                 double sum = 0.0;\r
 848                 double[] ev_array = o_ev.getArray();\r
 849                 int ev_clm = o_ev.getClm();\r
 850                 // For順変更禁止\r
 851                 for (i = 0; i < ev_clm; i++) {// for(int i = 0; i < ev->clm; i++ ){\r
 852                         sum += ev_array[i];// sum += ev->v[i];\r
 853                 }\r
 854                 // For順変更禁止\r
 855                 for (i = 0; i < ev_clm; i++) {// for(int i = 0; i < ev->clm; i++ ){\r
 856                         ev_array[i] /= sum;// ev->v[i] /= sum;\r
 857                 }\r
 858                 return;\r
 859         }\r
 860         \r
 861 \r
 862 \r
 863         private final NyARVec wk_vecTridiagonalize_vec = new NyARVec(0);\r
 864 \r
 865         private final NyARVec wk_vecTridiagonalize_vec2 = new NyARVec(0);\r
 866 \r
 867         /**\r
 868          * arVecTridiagonalize関数の代替品 a,d,e間で演算をしてる。何をどうしているかはさっぱりさっぱり\r
 869          * \r
 870          * @param a\r
 871          * @param d\r
 872          * @param e\r
 873          * @param i_e_start\r
 874          *            演算開始列(よくわからないけどarVecTridiagonalizeの呼び出し元でなんかしてる)\r
 875          * @return\r
 876          * @throws NyARException\r
 877          */\r
 878         private void vecTridiagonalize(NyARVec d, NyARVec e, int i_e_start)throws NyARException\r
 879         {\r
 880                 NyARVec vec = wk_vecTridiagonalize_vec;\r
 881                 // double[][] a_array=a.getArray();\r
 882                 double s, t, p, q;\r
 883                 int dim;\r
 884 \r
 885                 if (this.clm != this.row) {// if(a.getClm()!=a.getRow()){\r
 886                         throw new NyARException();\r
 887                 }\r
 888                 if (this.clm != d.getClm()) {// if(a.getClm() != d.clm){\r
 889                         throw new NyARException();\r
 890                 }\r
 891                 if (this.clm != e.getClm()) {// if(a.getClm() != e.clm){\r
 892                         throw new NyARException();\r
 893                 }\r
 894                 dim = this.getClm();\r
 895 \r
 896                 double[] d_vec, e_vec;\r
 897                 d_vec = d.getArray();\r
 898                 e_vec = e.getArray();\r
 899                 double[] a_vec_k;\r
 900 \r
 901                 for (int k = 0; k < dim - 2; k++) {\r
 902 \r
 903                         a_vec_k = this._m[k];\r
 904                         vec.setNewArray(a_vec_k, clm);// vec=this.getRowVec(k);//double[]\r
 905                                                                                         // vec_array=vec.getArray();\r
 906                         NyARException.trap("未チェックパス");\r
 907                         d_vec[k] = a_vec_k[k];// d.v[k]=vec.v[k];//d.set(k,v.get(k));\r
 908                                                                         // //d->v[k] = v[k];\r
 909 \r
 910                         // wv1.clm = dim-k-1;\r
 911                         // wv1.v = &(v[k+1]);\r
 912                         NyARException.trap("未チェックパス");\r
 913                         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
 914                         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
 915                                 continue;\r
 916                         }\r
 917 \r
 918                         for (int i = k + 1; i < dim; i++) {\r
 919                                 s = 0.0;\r
 920                                 for (int j = k + 1; j < i; j++) {\r
 921                                         NyARException.trap("未チェックのパス");\r
 922                                         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
 923                                 }\r
 924                                 for (int j = i; j < dim; j++) {\r
 925                                         NyARException.trap("未チェックのパス");\r
 926                                         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
 927                                 }\r
 928                                 NyARException.trap("未チェックのパス");\r
 929                                 d_vec[i] = s;// d.v[i]=s;//d->v[i] = s;\r
 930                         }\r
 931 \r
 932                         // wv1.clm = wv2.clm = dim-k-1;\r
 933                         // wv1.v = &(v[k+1]);\r
 934                         // wv2.v = &(d->v[k+1]);\r
 935                         a_vec_k = this._m[k];\r
 936                         vec.setNewArray(a_vec_k, clm);// vec=this.getRowVec(k);\r
 937                         // vec_array=vec.getArray();\r
 938                         NyARException.trap("未チェックパス");\r
 939                         t = vec.vecInnerproduct(d, k + 1) / 2;\r
 940                         for (int i = dim - 1; i > k; i--) {\r
 941                                 NyARException.trap("未チェックパス");\r
 942                                 p = a_vec_k[i];// p = v.get(i);//p = v[i];\r
 943                                 d_vec[i] -= t * p;\r
 944                                 q = d_vec[i];// d.v[i]-=t*p;q=d.v[i];//q = d->v[i] -= t*p；\r
 945                                 for (int j = i; j < dim; j++) {\r
 946                                         NyARException.trap("未チェックパス");\r
 947                                         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
 948                                 }\r
 949                         }\r
 950                 }\r
 951 \r
 952                 if (dim >= 2) {\r
 953                         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
 954                         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
 955                 }\r
 956 \r
 957                 if (dim >= 1) {\r
 958                         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
 959                 }\r
 960                 NyARVec vec2 = this.wk_vecTridiagonalize_vec2;\r
 961                 for (int k = dim - 1; k >= 0; k--) {\r
 962                         a_vec_k = this._m[k];\r
 963                         vec.setNewArray(a_vec_k, clm);// vec=this.getRowVec(k);//v =a.getPointer(k*dim);//v = &(a->m[k*dim]);\r
 964                         if (k < dim - 2) {\r
 965                                 for (int i = k + 1; i < dim; i++) {\r
 966                                         // wv1.clm = wv2.clm = dim-k-1;\r
 967                                         // wv1.v = &(v[k+1]);\r
 968                                         // wv2.v = &(a->m[i*dim+k+1]);\r
 969                                         vec2.setNewArray(this._m[i], clm);// vec2=this.getRowVec(i);\r
 970 \r
 971                                         t = vec.vecInnerproduct(vec2, k + 1);\r
 972                                         for (int j = k + 1; j < dim; j++) {\r
 973                                                 NyARException.trap("未チェックパス");\r
 974                                                 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
 975                                         }\r
 976                                 }\r
 977                         }\r
 978                         for (int i = 0; i < dim; i++) {\r
 979                                 a_vec_k[i] = 0.0;// v.set(i,0.0);//v[i] = 0.0;\r
 980                         }\r
 981                         a_vec_k[k] = 1;// v.set(k,1);//v[k] = 1;\r
 982                 }\r
 983                 return;\r
 984         }\r
 985 }