src/jp/nyatla/nyartoolkit/core/transmat/NyARTransMat.java

   1 /* \r
   2  * PROJECT: NyARToolkit (Extension)\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 edition ARToolKit class library.\r
  11  * Copyright (C)2008-2009 Ryo Iizuka\r
  12  *\r
  13  * This program is free software: you can redistribute it and/or modify\r
  14  * it under the terms of the GNU General Public License as published by\r
  15  * the Free Software Foundation, either version 3 of the License, or\r
  16  * (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 program.  If not, see <http://www.gnu.org/licenses/>.\r
  25  * \r
  26  * For further information please contact.\r
  27  *      http://nyatla.jp/nyatoolkit/\r
  28  *      <airmail(at)ebony.plala.or.jp> or <nyatla(at)nyatla.jp>\r
  29  * \r
  30  */\r
  31 package jp.nyatla.nyartoolkit.core.transmat;\r
  32 \r
  33 import jp.nyatla.nyartoolkit.NyARException;\r
  34 import jp.nyatla.nyartoolkit.core.param.*;\r
  35 import jp.nyatla.nyartoolkit.core.squaredetect.NyARSquare;\r
  36 import jp.nyatla.nyartoolkit.core.transmat.solver.*;\r
  37 import jp.nyatla.nyartoolkit.core.transmat.optimize.*;\r
  38 import jp.nyatla.nyartoolkit.core.transmat.rotmatrix.*;\r
  39 import jp.nyatla.nyartoolkit.core.types.*;\r
  40 import jp.nyatla.nyartoolkit.core.types.matrix.*;\r
  41 \r
  42 /**\r
  43  * This class calculates ARMatrix from square information and holds it. --\r
  44  * 変換行列を計算して、結果を保持するクラス。\r
  45  * \r
  46  */\r
  47 public class NyARTransMat implements INyARTransMat\r
  48 {       \r
  49         private final NyARDoublePoint2d _center=new NyARDoublePoint2d(0,0);\r
  50         private final NyARTransOffset _offset=new NyARTransOffset();\r
  51         private NyARPerspectiveProjectionMatrix _projection_mat_ref;\r
  52         protected NyARRotMatrix _rotmatrix;\r
  53         protected INyARTransportVectorSolver _transsolver;\r
  54         protected NyARPartialDifferentiationOptimize _mat_optimize;\r
  55 \r
  56 \r
  57         private NyARCameraDistortionFactor _ref_dist_factor;\r
  58 \r
  59         /**\r
  60          * 派生クラスで自分でメンバオブジェクトを指定したい場合はこちらを使う。\r
  61          *\r
  62          */\r
  63         protected NyARTransMat()\r
  64         {\r
  65                 //_calculator,_rotmatrix,_mat_optimizeをコンストラクタの終了後に\r
  66                 //作成して割り当ててください。\r
  67                 return;\r
  68         }\r
  69         public NyARTransMat(NyARParam i_param) throws NyARException\r
  70         {\r
  71                 final NyARCameraDistortionFactor dist=i_param.getDistortionFactor();\r
  72                 final NyARPerspectiveProjectionMatrix pmat=i_param.getPerspectiveProjectionMatrix();\r
  73                 this._transsolver=new NyARTransportVectorSolver(pmat,4);\r
  74                 //互換性が重要な時は、NyARRotMatrix_ARToolKitを使うこと。\r
  75                 //理屈はNyARRotMatrix_NyARToolKitもNyARRotMatrix_ARToolKitも同じだけど、少しだけ値がずれる。\r
  76                 this._rotmatrix = new NyARRotMatrix(pmat);\r
  77                 this._mat_optimize=new NyARPartialDifferentiationOptimize(pmat);\r
  78                 this._ref_dist_factor=dist;\r
  79                 this._projection_mat_ref=pmat;\r
  80         }\r
  81 \r
  82         public void setCenter(double i_x, double i_y)\r
  83         {\r
  84                 this._center.x= i_x;\r
  85                 this._center.y= i_y;\r
  86         }\r
  87 \r
  88 \r
  89 \r
  90 \r
  91         /**\r
  92          * 頂点順序をi_directionに対応して並べ替えます。\r
  93          * @param i_square\r
  94          * @param i_direction\r
  95          * @param o_sqvertex_ref\r
  96          * @param o_liner_ref\r
  97          */\r
  98         private final void initVertexOrder(NyARSquare i_square, int i_direction, NyARDoublePoint2d[] o_sqvertex_ref, NyARLinear[] o_liner_ref)\r
  99         {\r
 100                 //頂点順序を考慮した矩形の頂点情報\r
 101                 o_sqvertex_ref[0]= i_square.sqvertex[(4 - i_direction) % 4];\r
 102                 o_sqvertex_ref[1]= i_square.sqvertex[(5 - i_direction) % 4];\r
 103                 o_sqvertex_ref[2]= i_square.sqvertex[(6 - i_direction) % 4];\r
 104                 o_sqvertex_ref[3]= i_square.sqvertex[(7 - i_direction) % 4];    \r
 105                 o_liner_ref[0]=i_square.line[(4 - i_direction) % 4];\r
 106                 o_liner_ref[1]=i_square.line[(5 - i_direction) % 4];\r
 107                 o_liner_ref[2]=i_square.line[(6 - i_direction) % 4];\r
 108                 o_liner_ref[3]=i_square.line[(7 - i_direction) % 4];\r
 109                 return;\r
 110         }\r
 111 \r
 112 \r
 113         private final NyARDoublePoint2d[] __transMat_sqvertex_ref = new NyARDoublePoint2d[4];\r
 114         private final NyARDoublePoint2d[] __transMat_vertex_2d = NyARDoublePoint2d.createArray(4);\r
 115         private final NyARDoublePoint3d[] __transMat_vertex_3d = NyARDoublePoint3d.createArray(4);\r
 116         private final NyARLinear[] __transMat_linear_ref=new NyARLinear[4];\r
 117         private final NyARDoublePoint3d __transMat_trans=new NyARDoublePoint3d();\r
 118         \r
 119         /**\r
 120          * 頂点情報を元に、エラー閾値を計算します。\r
 121          * @param i_vertex\r
 122          */\r
 123         private double makeErrThreshold(NyARDoublePoint2d[] i_vertex)\r
 124         {\r
 125                 double a,b,l1,l2;\r
 126                 a=i_vertex[0].x-i_vertex[2].x;\r
 127                 b=i_vertex[0].y-i_vertex[2].y;\r
 128                 l1=a*a+b*b;\r
 129                 a=i_vertex[1].x-i_vertex[3].x;\r
 130                 b=i_vertex[1].y-i_vertex[3].y;\r
 131                 l2=a*a+b*b;\r
 132                 return (Math.sqrt(l1>l2?l1:l2))/200;\r
 133         }\r
 134         \r
 135         /**\r
 136          * double arGetTransMat( ARMarkerInfo *marker_info,double center[2], double width, double conv[3][4] )\r
 137          * \r
 138          * @param i_square\r
 139          * 計算対象のNyARSquareオブジェクト\r
 140          * @param i_direction\r
 141          * @param i_width\r
 142          * @return\r
 143          * @throws NyARException\r
 144          */\r
 145         public void transMat(final NyARSquare i_square, int i_direction, double i_width, NyARTransMatResult o_result_conv) throws NyARException\r
 146         {\r
 147                 final NyARDoublePoint2d[] sqvertex_ref = __transMat_sqvertex_ref;\r
 148                 final NyARLinear[] linear_ref=__transMat_linear_ref;\r
 149                 final NyARDoublePoint3d trans=this.__transMat_trans;\r
 150                 \r
 151                 double err_threshold=makeErrThreshold(i_square.sqvertex);\r
 152                 //計算用に頂点情報を初期化（順番調整）\r
 153                 initVertexOrder(i_square, i_direction, sqvertex_ref,linear_ref);\r
 154                 \r
 155                 //平行移動量計算機に、2D座標系をセット\r
 156                 NyARDoublePoint2d[] vertex_2d=this.__transMat_vertex_2d;\r
 157                 NyARDoublePoint3d[] vertex_3d=this.__transMat_vertex_3d;\r
 158                 this._ref_dist_factor.ideal2ObservBatch(sqvertex_ref, vertex_2d,4);             \r
 159                 this._transsolver.set2dVertex(vertex_2d,4);\r
 160                 \r
 161                 //基準矩形の3D座標系を作成\r
 162                 this._offset.setSquare(i_width,this._center);\r
 163 \r
 164                 //回転行列を計算\r
 165                 this._rotmatrix.initRotBySquare(linear_ref,sqvertex_ref);\r
 166                 \r
 167                 //回転後の3D座標系から、平行移動量を計算\r
 168                 this._rotmatrix.getPoint3dBatch(this._offset.vertex,vertex_3d,4);\r
 169                 this._transsolver.solveTransportVector(vertex_3d,trans);\r
 170                 \r
 171                 //計算結果の最適化(平行移動量と回転行列の最適化)\r
 172                 o_result_conv.error=this.optimize(this._rotmatrix, trans, this._transsolver,this._offset.vertex, vertex_2d,err_threshold);\r
 173                 \r
 174                 // マトリクスの保存\r
 175                 this.updateMatrixValue(this._rotmatrix, this._offset.point, trans,o_result_conv);\r
 176                 return;\r
 177         }\r
 178 \r
 179         /*\r
 180          * (non-Javadoc)\r
 181          * @see jp.nyatla.nyartoolkit.core.transmat.INyARTransMat#transMatContinue(jp.nyatla.nyartoolkit.core.NyARSquare, int, double, jp.nyatla.nyartoolkit.core.transmat.NyARTransMatResult)\r
 182          */\r
 183         public void transMatContinue(NyARSquare i_square, int i_direction, double i_width, NyARTransMatResult o_result_conv) throws NyARException\r
 184         {\r
 185                 final NyARDoublePoint2d[] sqvertex_ref = __transMat_sqvertex_ref;\r
 186                 final NyARLinear[] linear_ref=__transMat_linear_ref;\r
 187                 final NyARDoublePoint3d trans=this.__transMat_trans;\r
 188 \r
 189                 // io_result_convが初期値なら、transMatで計算する。\r
 190                 if (!o_result_conv.has_value) {\r
 191                         this.transMat(i_square, i_direction, i_width, o_result_conv);\r
 192                         return;\r
 193                 }\r
 194                 \r
 195                 //最適化計算の閾値を決定\r
 196                 double err_threshold=makeErrThreshold(i_square.sqvertex);\r
 197                 //計算用に頂点情報を初期化（順番調整）\r
 198                 initVertexOrder(i_square, i_direction, sqvertex_ref,linear_ref);\r
 199 \r
 200                 \r
 201                 //平行移動量計算機に、2D座標系をセット\r
 202                 NyARDoublePoint2d[] vertex_2d=this.__transMat_vertex_2d;\r
 203                 NyARDoublePoint3d[] vertex_3d=this.__transMat_vertex_3d;\r
 204                 this._ref_dist_factor.ideal2ObservBatch(sqvertex_ref, vertex_2d,4);             \r
 205                 this._transsolver.set2dVertex(vertex_2d,4);\r
 206                 \r
 207                 //基準矩形の3D座標系を作成\r
 208                 this._offset.setSquare(i_width,this._center);\r
 209 \r
 210                 //回転行列を計算\r
 211                 this._rotmatrix.initRotByPrevResult(o_result_conv);\r
 212                 \r
 213                 //回転後の3D座標系から、平行移動量を計算\r
 214                 this._rotmatrix.getPoint3dBatch(this._offset.vertex,vertex_3d,4);\r
 215                 this._transsolver.solveTransportVector(vertex_3d,trans);\r
 216 \r
 217                 //現在のエラーレートを計算しておく\r
 218                 double min_err=errRate(this._rotmatrix,trans, this._offset.vertex, vertex_2d,4,vertex_3d);\r
 219                 NyARDoubleMatrix33 rot=this.__rot;\r
 220                 //エラーレートが前回のエラー値より閾値分大きかったらアゲイン\r
 221                 if(min_err<o_result_conv.error+err_threshold){\r
 222                         rot.setValue(this._rotmatrix);\r
 223                         //最適化してみる。\r
 224                         for (int i = 0;i<5; i++) {\r
 225                                 //変換行列の最適化\r
 226                                 this._mat_optimize.modifyMatrix(rot, trans, this._offset.vertex, vertex_2d, 4);\r
 227                                 double err=errRate(rot,trans,this._offset.vertex, vertex_2d,4,vertex_3d);\r
 228                                 //System.out.println("E:"+err);\r
 229                                 if(min_err-err<err_threshold/2){\r
 230                                         //System.out.println("BREAK");\r
 231                                         break;\r
 232                                 }\r
 233                                 this._transsolver.solveTransportVector(vertex_3d, trans);\r
 234                                 this._rotmatrix.setValue(rot);\r
 235                                 min_err=err;\r
 236                         }\r
 237                         this.updateMatrixValue(this._rotmatrix, this._offset.point, trans,o_result_conv);\r
 238                 }else{\r
 239                         //回転行列を計算\r
 240                         this._rotmatrix.initRotBySquare(linear_ref,sqvertex_ref);\r
 241                         \r
 242                         //回転後の3D座標系から、平行移動量を計算\r
 243                         this._rotmatrix.getPoint3dBatch(this._offset.vertex,vertex_3d,4);\r
 244                         this._transsolver.solveTransportVector(vertex_3d,trans);\r
 245                         \r
 246                         //計算結果の最適化(平行移動量と回転行列の最適化)\r
 247                         min_err=this.optimize(this._rotmatrix, trans, this._transsolver,this._offset.vertex, vertex_2d,err_threshold);\r
 248                         this.updateMatrixValue(this._rotmatrix, this._offset.point, trans,o_result_conv);\r
 249                 }\r
 250                 o_result_conv.error=min_err;\r
 251                 return;\r
 252         }\r
 253         private NyARDoubleMatrix33 __rot=new NyARDoubleMatrix33();\r
 254         private double optimize(NyARRotMatrix io_rotmat,NyARDoublePoint3d io_transvec,INyARTransportVectorSolver i_solver,NyARDoublePoint3d[] i_offset_3d,NyARDoublePoint2d[] i_2d_vertex,double i_err_threshold) throws NyARException\r
 255         {\r
 256                 //System.out.println("START");\r
 257                 NyARDoublePoint3d[] vertex_3d=this.__transMat_vertex_3d;\r
 258                 //初期のエラー値を計算\r
 259                 double min_err=errRate(io_rotmat, io_transvec, i_offset_3d, i_2d_vertex,4,vertex_3d);\r
 260                 NyARDoubleMatrix33 rot=this.__rot;\r
 261                 rot.setValue(io_rotmat);\r
 262                 for (int i = 0;i<5; i++) {\r
 263                         //変換行列の最適化\r
 264                         this._mat_optimize.modifyMatrix(rot, io_transvec, i_offset_3d, i_2d_vertex, 4);\r
 265                         double err=errRate(rot,io_transvec, i_offset_3d, i_2d_vertex,4,vertex_3d);\r
 266                         //System.out.println("E:"+err);\r
 267                         if(min_err-err<i_err_threshold){\r
 268                                 //System.out.println("BREAK");\r
 269                                 break;\r
 270                         }\r
 271                         i_solver.solveTransportVector(vertex_3d, io_transvec);\r
 272                         io_rotmat.setValue(rot);\r
 273                         min_err=err;\r
 274                 }\r
 275                 //System.out.println("END");\r
 276                 return min_err;\r
 277         }\r
 278         \r
 279         //エラーレート計算機\r
 280         public double errRate(NyARDoubleMatrix33 io_rot,NyARDoublePoint3d i_trans, NyARDoublePoint3d[] i_vertex3d, NyARDoublePoint2d[] i_vertex2d,int i_number_of_vertex,NyARDoublePoint3d[] o_rot_vertex) throws NyARException\r
 281         {\r
 282                 NyARPerspectiveProjectionMatrix cp = this._projection_mat_ref;\r
 283                 final double cp00=cp.m00;\r
 284                 final double cp01=cp.m01;\r
 285                 final double cp02=cp.m02;\r
 286                 final double cp11=cp.m11;\r
 287                 final double cp12=cp.m12;\r
 288 \r
 289                 double err=0;\r
 290                 for(int i=0;i<i_number_of_vertex;i++){\r
 291                         double x3d,y3d,z3d;\r
 292                         o_rot_vertex[i].x=x3d=io_rot.m00*i_vertex3d[i].x+io_rot.m01*i_vertex3d[i].y+io_rot.m02*i_vertex3d[i].z;\r
 293                         o_rot_vertex[i].y=y3d=io_rot.m10*i_vertex3d[i].x+io_rot.m11*i_vertex3d[i].y+io_rot.m12*i_vertex3d[i].z;\r
 294                         o_rot_vertex[i].z=z3d=io_rot.m20*i_vertex3d[i].x+io_rot.m21*i_vertex3d[i].y+io_rot.m22*i_vertex3d[i].z;\r
 295                         x3d+=i_trans.x;\r
 296                         y3d+=i_trans.y;\r
 297                         z3d+=i_trans.z;\r
 298                         \r
 299                         //射影変換\r
 300                         double x2d=x3d*cp00+y3d*cp01+z3d*cp02;\r
 301                         double y2d=y3d*cp11+z3d*cp12;\r
 302                         double h2d=z3d;\r
 303                         \r
 304                         //エラーレート計算\r
 305                         double t1=i_vertex2d[i].x-x2d/h2d;\r
 306                         double t2=i_vertex2d[i].y-y2d/h2d;\r
 307                         err+=t1*t1+t2*t2;\r
 308                         \r
 309                 }\r
 310                 return err/i_number_of_vertex;\r
 311         }               \r
 312         \r
 313         \r
 314         \r
 315         /**\r
 316          * パラメータで変換行列を更新します。\r
 317          * \r
 318          * @param i_rot\r
 319          * @param i_off\r
 320          * @param i_trans\r
 321          */\r
 322         public void updateMatrixValue(NyARRotMatrix i_rot, NyARDoublePoint3d i_off, NyARDoublePoint3d i_trans,NyARTransMatResult o_result)\r
 323         {\r
 324                 o_result.m00=i_rot.m00;\r
 325                 o_result.m01=i_rot.m01;\r
 326                 o_result.m02=i_rot.m02;\r
 327                 o_result.m03=i_rot.m00 * i_off.x + i_rot.m01 * i_off.y + i_rot.m02 * i_off.z + i_trans.x;\r
 328 \r
 329                 o_result.m10 = i_rot.m10;\r
 330                 o_result.m11 = i_rot.m11;\r
 331                 o_result.m12 = i_rot.m12;\r
 332                 o_result.m13 = i_rot.m10 * i_off.x + i_rot.m11 * i_off.y + i_rot.m12 * i_off.z + i_trans.y;\r
 333 \r
 334                 o_result.m20 = i_rot.m20;\r
 335                 o_result.m21 = i_rot.m21;\r
 336                 o_result.m22 = i_rot.m22;\r
 337                 o_result.m23 = i_rot.m20 * i_off.x + i_rot.m21 * i_off.y + i_rot.m22 * i_off.z + i_trans.z;\r
 338 \r
 339                 o_result.has_value = true;\r
 340                 return;\r
 341         }       \r
 342 }\r